Release v1.2.0

.gitignore

@@ -2,6 +2,7 @@
 __pycache__/
 *.py[cod]
 *$py.class
+*.drawio.bkp
 
 # C extensions
 *.so

CHANGELOG.md

@@ -1,6 +1,23 @@
 # Changelog
 All new changes are documented here.
 
+## [v1.2.0]
+### Added
+- `BGC` workflow detecting regions with high gene conversion rate in the genome using the alignments directly
+- `NEUTRAL_MODEL` workflow producing the neutral model needed for both `CONSTRAINED` and `BGC`
+
+### Changed
+- Optional `cactus-hal2maf` by using `--cactus_hal2maf` (default use legacy hal2maf)
+- Updated phast to v1.9.9
+- `mutyper_variant` process is now more streamlined
+- `CONSTRAINED`, `BGC`, `NEUTRAL_MODEL` and `IBD` workflow can now run when `--vcf_is_filtered` is passed
+- K-mer counting from the ancestral genome is now done using [smakcr](https://github.com/julibeg/smakcr)
+- Faster parallel processing of intervals in `NEUTRAL_MODEL`, `CONSTRAINED` and `BGC` workflows
+
+### Fixed
+- IBD workflow now runs to completion with demo data
+- Workflow breaking when all SDM files for a breed are empty
+
 ## [v1.1.6a]
 ### Fixed
 - Broken configuration file

Dockerfile

@@ -15,13 +15,12 @@ RUN mamba install -n nspectra -y \
     bioconda::bedtools bioconda::plink=1.90
 RUN mamba install -n nspectra -y bioconda::ucsc-twobitinfo bioconda::ucsc-fatotwobit \
     bioconda::ucsc-wigtobigwig bioconda::ucsc-bigwigtobedgraph
-RUN mamba install -n nspectra -y bioconda::vcflib==1.0.3
 RUN mamba install -n nspectra -y bioconda::shapeit4
 RUN mamba install -n nspectra -y bioconda::shapeit5
 RUN mamba install -n nspectra -y \
     bioconda::perl-bioperl \
-    bioconda::phast=1.5 bioconda::mutyper bioconda::perl-bio-db-hts \
-    bioconda::tabix bioconda::tabixpp=1.1.0
+    bioconda::phast=1.9.9 bioconda::mutyper bioconda::perl-bio-db-hts \
+    bioconda::tabix bioconda::vcflib
 RUN mamba install -n nspectra -y \
     conda-forge::r-base>=4.1.0 \
     conda-forge::r-cowplot \
@@ -46,6 +45,14 @@ RUN conda-pack -n nspectra -o /tmp/env.tar && \
 RUN /venv/bin/conda-unpack
 
 
+# Compile and use Smakcr
+FROM rust:1.97.1-alpine AS cargo_build
+WORKDIR /opt/
+RUN apk add --no-cache git
+RUN git clone https://github.com/julibeg/smakcr
+RUN cd smakcr && cargo build --release
+
+
 # The runtime-stage image; we can use Debian as the
 # base image since the Conda env also includes Python
 # for us.
@@ -60,6 +67,7 @@ RUN chmod a+x /usr/local/bin/datasets
 
 # Copy /venv from the previous stage:
 COPY --from=build /venv /venv
+COPY --from=cargo_build /opt/smakcr/target/release/smakcr /usr/local/bin/smakcr
 
 # When image is run, run the code with the environment
 # activated:

README.md

@@ -63,6 +63,21 @@ Some steps of the workflow are very computationally intensive, and can benefit f
 ## Input pre-filtering
 We recommend to pre-filter the vcf file to obtain samples with a reasonable coverage (>8-10 mean DP) and with variants with low call rate removed from the dataset (CCR > 90%) and minor allele count MAC >= 2.
 
+## Remove constrained regions
+> NOTE: this is an experimental feature, and needs further testing
+It is possible to exclude conserved elements identified from the whole genome alignments directly. To do so, run the workflow with `--remove_constrained`. To run this feature, you need to provide a BED file with the exon with `--exon_bed`.
+
+## Biased Gene Conversion (gBGC) removal
+> NOTE: this is an experimental feature, and needs further testing
+It is possible to exclude gBGC regions identified from the whole genome alignments using [phastBias](http://compgen.cshl.edu/phast/help-pages/phastBias.txt). To do so, run the workflow with `--bgc`. To run this feature, you need to provide a BED file with the exon with `--exon_bed`.
+
+### Create an exon bed file
+You can create an exon bed file from your own GFF3 annotation using `gff3ToGenePred` followed by `genePredToBed` (see [here](https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/)):
+```
+gff3ToGenePred annotation.gff3 annotation.gp
+genePredToBed annotation.gp exons.bed
+```
+
 # How to cite
 To cite nSPECTRa refer to:
 

conf/base.config

@@ -35,12 +35,12 @@ process {
     time = { 12.h * task.attempt }
   }
   withLabel: medium_vlargemem{
-    cpus = 4
+    cpus = 8
     memory = { 512.GB * task.attempt }
-    time = { 12.h * task.attempt }
+    time = { 23.h * task.attempt }
   }
   withLabel: medium_multi{
-    cpus = 4
+    cpus = 8
     memory = { 64.GB * task.attempt }
     time = { 12.h * task.attempt }
   }

conf/demo.config

@@ -4,7 +4,8 @@ params {
     idx = "${projectDir}/test/demo.vcf.gz.tbi"
     hal = "${projectDir}/test/demo.hal"
     gff = "${projectDir}/test/chr24.gff.gz"
-    reference = 'Bbu.GCF_003121395.1'
+    exon_bed = "${projectDir}/test/exons.bed"
+    reference = 'Bbu_GCF_003121395'
     target = "Inner2"
     outdir = 'results'
     species = 'bubalus_bubalis'
@@ -20,6 +21,8 @@ params {
     mutyper = true
     sdm = true
     relate = false
+    bgc = true
+    remove_constrained = true
     // Preprocessing
     vep = true
     vcf_is_filtered = false
@@ -31,14 +34,14 @@ params {
     gone = false
     refinedibd = false
     mergeibd = false
-    compute_ibd = false
+    compute_ibd = true
 
     // Settings
     k = '3,5,7'
     ref_min_size = false
 
     // SDM
-    chunk_size = 100000
+    chunk_size = 10000000
 
     // Effective pop size
     neval = false

include/process/hal2maf.nf

@@ -17,6 +17,26 @@ process hal2maf {
 
 
     script:
+    if (params.cactus_hal2maf)
+    """
+    export HOME=\$PWD
+    mkdir toil-work
+    mkdir toil-coord
+    toil config test.yaml
+    cactus-hal2maf \
+        ./js \
+        ${HAL} \
+        ${params.reference}_${params.target}.maf \
+        --config test.yaml \
+        --refGenome ${params.reference} \
+        --targetGenomes ${params.target} \
+        --batchCores 1 \
+        --defaultCores 1 \
+        --workDir ./toil-work \
+        --coordinationDir ./toil-coord \
+        --cleanWorkDir onSuccess
+    """
+    else
     """
     hal2maf \
         --refGenome ${params.reference} \

include/process/ibd.nf

@@ -4,18 +4,21 @@ process ibd {
     label "medium"
 
     input:
-    path vcf
-    path tbi
+    tuple val(contig), path(vcf), path(tbi), path(map)
     path refinedibd
-    val contig
 
     output:
-    tuple val(contig), path("ibd.${contig}.ibd.gz")
+    tuple val(contig), path(vcf), path(tbi), path(map), path("IBD.${contig}.ibd.gz"), emit: ibd
+    tuple val(contig), path(vcf), path(tbi), path(map), path("IBD.${contig}.hbd.gz"), emit: hbd
 
     script:
     """
     javamem=`python -c "import sys; maxmem=int(sys.argv[1]); print( maxmem - int(maxmem * .1) )" ${ task.memory.toGiga() }`
-    java -Xmx\${javamem}G -jar ${refinedibd} ${params.refined_ibd_params} \
+    java -Xmx\${javamem}G -jar ${refinedibd} \
+            window=${params.refined_ibd_window} \
+            trim=${params.refined_ibd_trim} \
+            length=${params.refined_ibd_length} \
+            lod=${params.refined_ibd_lod} \
             gt=${vcf} \
             out=IBD.${contig} \
             chrom=${contig} \
@@ -31,15 +34,15 @@ process ibd {
 process make_map {
 
     input:
-    path vcf
-    path tbi
+    tuple val(chrom), path(vcf), path(tbi)
 
     output:
-    path "${vcf.simpleName}.map"
+
+    tuple val(chrom), path(vcf), path(tbi), path("${vcf.simpleName}.map")
 
     script:
     """
-    bcftools query -f '%CHROM\t%ID\t%POS\n' ${vcf} | awk 'BEGIN{OFS="\t"}; {print \$1,\$2,\$3/1000000,\$3}'> ${vcf.simpleName}.map
+    bcftools query -f '%CHROM\t%ID\t%POS\n' ${vcf} | awk 'BEGIN{OFS="\t"}; {print \$1,\$2,\$3 * ${params.recombination_rate},\$3}'> ${vcf.simpleName}.map
     """
 
     stub:
@@ -51,30 +54,53 @@ process make_map {
 process merge_ibd {
     tag "ibd ${contig}"
     label "medium"
+    publishDir "${params.outdir}/IBD/single", mode: "${params.publish_dir_mode}", overwrite: true
 
     input:
-    tuple val(contig), path(ibd)
+    tuple val(contig), path(vcf), path(tbi), path(map), path(ibd)
     path merge_ibd
-    path vcf
-    path tbi
-    path map
+    val output_root
 
     output:
-    path "ibd.${contig}.ibd.gz"
+    path "${output_root}.${contig}.merged.ibd.gz"
 
     script:
     """
     javamem=`python -c "import sys; maxmem=int(sys.argv[1]); print( maxmem - int(maxmem * .1) )" ${ task.memory.toGiga() }`
-    java -Xmx\${javamem}G -jar ${merge_ibd} \
-            ${ibd} \
+    zcat ${ibd} | java -Xmx\${javamem}G -jar ${merge_ibd} \
             ${vcf} \
             ${map} \
-            ${params.merge_ibd_params} \
-            IBD.${contig}.merge
+            ${params.merge_ibd_gap} ${params.merge_ibd_discordance} | \
+            bgzip -c > ${output_root}.${contig}.merged.ibd.gz
     """
 
     stub:
     """
-    touch ibd.${contig}.ibd.gz
+    touch ${output_root}.${contig}.merged.ibd.gz
     """
-}
+}
+
+process zcatsort {
+    publishDir { "${params.outdir}/${OUTDIR}" }, mode: "${params.publish_dir_mode}", overwrite: true
+    label "small"
+
+    input:
+    path INFILES
+    val OUTNAME
+    val OUTDIR
+    val SORTEXPR
+
+    output:
+    path "${OUTNAME}"
+
+
+    script:
+    """
+    zcat ${INFILES} | sort ${SORTEXPR} | bgzip -c > ${OUTNAME}
+    """
+
+    stub:
+    """
+    touch ${OUTNAME}
+    """
+}

include/process/mutyper.nf

@@ -31,20 +31,16 @@ process mutyper_variant {
     if (params.annotation)
     """
     echo "Run mutyper (variants)"
-    bcftools view --threads ${task.cpus} -v snps -r ${chrom} -m2 -M2 ${vcf} | \
-        bedtools intersect -header -v -a - -b ${masks_ch} | \
-        sed 's/_pilon//g' | \
-        vcffixup - ${vcftools_filter} |\
-        mutyper variants --k ${k} --strand_file ${params.annotation} ${ancfasta} - | \
-        bgzip -c > mutationSpectra_${params.species.capitalize()}_${chrom}_${k}.vcf.gz &&
+    bcftools view --threads ${task.cpus} -v snps -r ${chrom} -T ^${masks_ch} -m2 -M2 ${vcf} | \\
+        vcffixup - ${vcftools_filter} |\\
+        mutyper variants --k ${k} --strand_file ${params.annotation} ${ancfasta} - | \\
+        bgzip -c > mutationSpectra_${params.species.capitalize()}_${chrom}_${k}.vcf.gz && \\
         tabix -p vcf mutationSpectra_${params.species.capitalize()}_${chrom}_${k}.vcf.gz
     """
     else
     """
     echo "Run mutyper (variants)"
-    bcftools view --threads ${task.cpus} -v snps -r ${chrom} -m2 -M2 ${vcf} | \
-        bedtools intersect -header -v -a - -b ${masks_ch} | \
-        sed 's/_pilon//g' | \
+    bcftools view --threads ${task.cpus} -v snps -r ${chrom} -T ^${masks_ch} -m2 -M2 ${vcf} | \
         vcffixup - ${vcftools_filter} |\
         mutyper variants --k ${k} ${ancfasta} - | \
         bgzip -c > mutationSpectra_${params.species.capitalize()}_${chrom}_${k}.vcf.gz &&
@@ -317,9 +313,14 @@ process ksfs {
     """
     echo "Run mutyper (ksfs)"
     bcftools query -l ${vcf} > all_samples.txt
-    grep -f all_samples.txt ${samplelist} > keep.txt
-    bcftools view --threads ${task.cpus} -S keep.txt ${vcf} | \
-        mutyper ksfs - > ksfs_${samplename}_${k}.tsv
+    if ! grep -f all_samples.txt ${samplelist}; then
+        # Create an empty Ksfs file for empty populations.
+        touch ksfs_${samplename}_${k}.tsv
+    else
+        grep -f all_samples.txt ${samplelist} > keep.txt    
+        bcftools view --threads ${task.cpus} -S keep.txt ${vcf} | \
+            mutyper ksfs - > ksfs_${samplename}_${k}.tsv
+    fi
     """
 
     stub:
@@ -343,7 +344,7 @@ process kmercount {
 
     script:
     """
-    kmer_count ${ancfa} ${k} > ${params.species.capitalize()}.K${k}.txt
+    smakcr -k ${k} -o ${params.species.capitalize()}.K${k}.txt ${ancfa}
     """
 
     stub: