diff --git a/.gitignore b/.gitignore index b6e4761..e3e7acd 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ __pycache__/ *.py[cod] *$py.class +*.drawio.bkp # C extensions *.so diff --git a/CHANGELOG.md b/CHANGELOG.md index fe823ff..dbf9079 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Dockerfile b/Dockerfile index ebe6f14..f645517 100755 --- a/Dockerfile +++ b/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: diff --git a/README.md b/README.md index 2511755..bb024bc 100755 --- a/README.md +++ b/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: diff --git a/conf/base.config b/conf/base.config index 60ee4b0..f0f5840 100755 --- a/conf/base.config +++ b/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 } } diff --git a/conf/demo.config b/conf/demo.config index 8d5bb42..c0ffc68 100755 --- a/conf/demo.config +++ b/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 diff --git a/include/process/hal2maf.nf b/include/process/hal2maf.nf index 7cdd18f..4183e85 100755 --- a/include/process/hal2maf.nf +++ b/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} \ diff --git a/include/process/ibd.nf b/include/process/ibd.nf index 7db6660..29ad88f 100755 --- a/include/process/ibd.nf +++ b/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 """ -} \ No newline at end of file +} + +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} + """ +} diff --git a/include/process/mutyper.nf b/include/process/mutyper.nf index e4a0cdc..53621d5 100755 --- a/include/process/mutyper.nf +++ b/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: diff --git a/include/process/phylop.nf b/include/process/phylop.nf index 1a695b5..cf2b13d 100755 --- a/include/process/phylop.nf +++ b/include/process/phylop.nf @@ -30,6 +30,56 @@ process hal4d { """ } +process split_4d_bed { + tag "h4d" + label "small" + + input: + tuple path(hal), val(GENOMES), path(neutral_bed) + + output: + tuple path(hal), val(GENOMES), path("*.4d.bed") + + script: + """ + sort -k1,1 -k 2,2n ${neutral_bed} | \ + bedtools merge -i - | \ + awk -v chunk_sz=10000000 'BEGIN {CNT=0; ITV=1}; CNT "interval_"ITV".bed"; CNT+=\$3-\$2}; CNT>=chunk_sz {ITV+=1; print \$0 > "interval_"ITV".bed"; CNT=\$3-\$2}' + for internal_bed in interval_*.bed; do + bname=\$( basename \${internal_bed} .bed ) + bedtools intersect -u -a ${neutral_bed} -b \$internal_bed > \${bname}.4d.bed + done + """ + + stub: + """ + touch interval_4d_1.bed + touch interval_4d_2.bed + touch interval_4d_3.bed + """ +} + +process collect_maf { + tag "h4d" + label "small" + + input: + tuple val(GENOMES), path(mafs) + + output: + tuple val(GENOMES), path("4d.maf") + + script: + """ + awk 'NR>1 || NR==1' ${mafs} > 4d.maf + """ + + stub: + """ + touch 4d.maf + """ +} + process halTree { tag "pFit" publishDir "${params.outdir}/PHYLOP/NEUTRAL", mode: "${params.publish_dir_mode}", overwrite: true @@ -81,21 +131,42 @@ process phyloFit { process make4dmaf { tag "maf" publishDir "${params.outdir}/PHYLOP/MAF", mode: "${params.publish_dir_mode}", overwrite: true - label "medium_vlargemem" + label "medium_largemem" conda {params.enable_conda ? "${baseDir}/envs/phast_environment.yml" : null} container { params.cactus_version ? "quay.io/comparative-genomics-toolkit/cactus:${params.cactus_version}" : "quay.io/comparative-genomics-toolkit/cactus:latest" } input: - tuple path(hal), path(bedfile), val(GENOMES) + tuple path(hal), val(GENOMES), path(bedfile) output: - tuple path(hal), path("4d.maf"), val(GENOMES) + tuple path(hal), val(GENOMES), path("${bedfile.simpleName}.4d.maf") script: + if (params.cactus_hal2maf) + """ + export HOME=\$PWD + toil config test.yaml + mkdir toil-work + mkdir toil-coord + cactus-hal2maf \ + ./js \ + ${hal} \ + ${bedfile.simpleName}.4d.maf \ + --config test.yaml \ + --refGenome ${params.reference} \ + --bedRanges ${bedfile} \ + --targetGenomes "${GENOMES}" \ + --batchCores ${task.cpus} \ + --defaultCores ${task.cpus} \ + --workDir ./toil-work \ + --coordinationDir ./toil-coord \ + --cleanWorkDir onSuccess + """ + else """ - hal2mafMP.py \ + hal2maf \ ${hal} \ - 4d.maf \ + ${bedfile.simpleName}.4d.maf \ --noDupes \ --targetGenomes "${GENOMES}" \ --numProc ${task.cpus} \ @@ -142,7 +213,7 @@ process msa_view { conda {params.enable_conda ? "${baseDir}/envs/phast_environment.yml" : null} input: - tuple path(hal), path(maf), val(GENOMES) + tuple val(GENOMES), path(maf) output: path "${maf.simpleName}.ss" @@ -195,21 +266,20 @@ process phyloP { container { params.cactus_version ? "quay.io/comparative-genomics-toolkit/cactus:${params.cactus_version}" : "quay.io/comparative-genomics-toolkit/cactus:latest" } input: - tuple val(n), val(chr), path(hal), path(model) + tuple path(hal), path(BED) + path model output: - tuple val(n), val(chr), path(hal), path(model), path("phylop_${chr}.wig"), path("${params.reference}.sizes") + path "${BED.simpleName}.wig" script: """ - halPhyloPMP.py \ + halPhyloP \ ${hal} \ ${params.reference} \ ${model} \ - phylop_${chr}.wig \ - --refSequence ${chr} \ - --chromSizes ${params.reference}.sizes \ - --numProc ${task.cpus} + ${BED.simpleName}.wig \ + --refBed ${BED} """ stub: @@ -225,7 +295,7 @@ process wig2bedgraph { afterScript "rm ${wig.baseName}.bw" input: - tuple val(n), val(chr), path(hal), path(model), path(wig), path(sizes) + tuple path(wig), path(sizes) output: path "${wig.baseName}.bed" @@ -270,19 +340,21 @@ process bedtobigwig { process combine_bed { tag "bed" - publishDir "${params.outdir}/PHYLOP/BED", mode: "${params.publish_dir_mode}", overwrite: true + publishDir { "${params.outdir}/${outdir}/BED" }, mode: "${params.publish_dir_mode}", overwrite: true label "largemem" conda {params.enable_conda ? "${baseDir}/envs/phast_environment.yml" : null} input: path beds + val outname + val outdir output: - path "phylop.bed" + path "${outname}.bed" script: """ - cat ${beds} | bedtools sort -i - > phylop.bed + cat ${beds} | sort -k 1,1 -k2,2n --parallel ${task.cpus} - > ${outname}.bed """ stub: @@ -336,30 +408,294 @@ process extract_conserved { } -process vcf_drop_conserved { +process vcf_drop_intervals { tag "filt" - publishDir "${params.outdir}/PHYLOP/VCF", mode: "${params.publish_dir_mode}", overwrite: true + publishDir { "${params.outdir}/${OUTDIR}/VCF" }, mode: "${params.publish_dir_mode}", overwrite: true label "medium" input: - path vcf - path tbi + tuple val(chrom), path(vcf), path(tbi) path bed + val tag + val OUTDIR output: - path "${vcf.simpleName}.non-conserved.vcf.gz", emit: vcf - path "${vcf.simpleName}.non-conserved.vcf.gz.tbi", emit: tbi + tuple val(chrom), path("${vcf.simpleName}.${tag}.vcf.gz"), path("${vcf.simpleName}.${tag}.vcf.gz.tbi") script: """ - bedtools intersect -header -v -a ${vcf} -b ${bed} | bgzip -c > ${vcf.simpleName}.non-conserved.vcf.gz - tabix -p vcf ${vcf.simpleName}.non-conserved.vcf.gz + bedtools intersect -header -v -a ${vcf} -b ${bed} | \ + bgzip -@ ${task.cpus > 1 ? task.cpus - 1 : 1} -c > ${vcf.simpleName}.${tag}.vcf.gz + tabix -p vcf ${vcf.simpleName}.${tag}.vcf.gz + bcftools stats ${vcf.simpleName}.${tag}.vcf.gz > ${vcf.simpleName}.${tag}.bcftools_stats """ stub: """ - touch ${vcf.simpleName}.non-conserved.vcf.gz - touch ${vcf.simpleName}.non-conserved.vcf.gz.tbi + touch ${vcf.simpleName}.${tag}.vcf.gz + touch ${vcf.simpleName}.${tag}.vcf.gz.tbi + """ +} + +// Processes for the phastBias methods +process GENOME_INTERVALS { + conda "bioconda::pysam=0.22.1" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/pysam:0.22.1--py39hcada746_0' : + 'quay.io/biocontainers/pysam:0.22.1--py39hcada746_0' }" + + input: + path sizes + + output: + path "intervals_*.bed" + + script: + """ + #!/usr/bin/env python + import pysam + import re + + # Load the sizes and order them + sizes_fn = open("${sizes}") + sizes = [] + for line in sizes_fn: + line = line.strip().split() + sizes.append((int(line[1]), line[0])) + sizes.sort(reverse=True) + + # Now, we can create the chunks + n = 1 + target_size = ${params.chunk_size} + proc_size = 0 + tmp_list = [] + tot_len = 0 + for (seq_len, seq_id) in sizes: + if seq_len > target_size: + for i in range(0, seq_len, target_size): + with open(f"intervals_{n}.bed", "w") as bedfile: + start = i + end = min(i + target_size, seq_len) + bedfile.write(f'{seq_id}\\t{start}\\t{end}\\n') + n+=1 + else: + if tot_len >= target_size: + with open(f"intervals_{n}.bed", "w") as bedfile: + for line in tmp_list: + bedfile.write(line) + n+=1 + tmp_list = [] + tot_len = 0 + else: + tmp_list.append(f'{seq_id}\\t0\\t{seq_len}\\n') + tot_len += seq_len + if len(tmp_list) > 0: + with open(f"intervals_{n}.bed", "w") as bedfile: + for line in tmp_list: + bedfile.write(line) """ } +process create_maf { + tag "hal2maf" + publishDir "${params.outdir}/MAF", mode: "${params.publish_dir_mode}", overwrite: true + container { params.cactus_version ? "quay.io/comparative-genomics-toolkit/cactus:${params.cactus_version}" : "quay.io/comparative-genomics-toolkit/cactus:latest" } + + label "largemem" + + input: + tuple path(HAL), path(BED) + val GENOMES + + output: + tuple path("${BED.simpleName}.maf"), path(BED) + + + script: + if (params.cactus_hal2maf) + """ + export HOME=\$PWD + mkdir toil-work + mkdir toil-coord + toil config test.yaml + cactus-hal2maf \ + ./js \ + ${HAL} \ + ${BED.simpleName}.maf \ + --config test.yaml \ + --refGenome ${params.reference} \ + --noAncestors \ + --bedRanges ${BED} \ + --batchCores 1 \ + --outType single \ + --defaultCores 1 \ + --workDir ./toil-work \ + --coordinationDir ./toil-coord \ + --cleanWorkDir onSuccess + """ + else + """ + hal2maf \ + --noAncestors \ + --refGenome ${params.reference} \ + --targetGenomes "${GENOMES}" \ + --refTargets ${BED} \ + --hdf5InMemory ${HAL} alignments.${CHROM}.maf + """ + + stub: + """ + touch alignments.${CHROM}.maf + """ +} + + + +process phastBias { + tag "phastBias" + publishDir "${params.outdir}/PHAST/gBGC", mode: "${params.publish_dir_mode}", overwrite: true + label "largemem" + conda {params.enable_conda ? "${baseDir}/envs/phast_environment.yml" : null} + + input: + tuple path(maf), path(BED) + path model + + output: + path "${maf.simpleName}.wig", emit: wig + path "${maf.simpleName}.tracts.bed", emit: tracts_gff + path "${maf.simpleName}.informative.bed", emit: informative_gff + path "${maf.simpleName}.tracts.bed", emit: tracts_bed + path "${maf.simpleName}.informative.bed", emit: informative_bed + + script: + """ + CHROM=\$( head -1 $BED | cut -f1 ) + phastBias \ + --informative-fn ${maf.simpleName}.informative.gff \ + --output-tracts ${maf.simpleName}.tracts.gff \ + ${maf} \ + ${model} \ + ${params.reference} | \ + sed "s/chrom=${params.reference}/chrom=\$CHROM/" > ${maf.simpleName}.wig + awk -v var=\$CHROM '{OFS="\\t"; print var, \$4, \$5, \$3}' ${maf.simpleName}.informative.gff > ${maf.simpleName}.informative.bed + awk -v var=\$CHROM '{OFS="\\t"; print var, \$4, \$5, \$3}' ${maf.simpleName}.tracts.gff > ${maf.simpleName}.tracts.bed + """ + + stub: + """ + touch ${maf.simpleName}.wig + touch ${maf.simpleName}.gff + """ +} + +process halSize { + tag "medium_mem" + publishDir "${params.outdir}/MAF", mode: "${params.publish_dir_mode}", overwrite: true + container { params.cactus_version ? "quay.io/comparative-genomics-toolkit/cactus:${params.cactus_version}" : "quay.io/comparative-genomics-toolkit/cactus:latest" } + + label "largemem" + + input: + path HAL + + output: + path "${params.reference}.sizes" + + + script: + """ + halStats --sequenceStats ${params.reference} ${HAL} | \ + awk 'BEGIN{FS=","}; NR>1 && \$1!="" {print \$1"\\t"\$2}' > ${params.reference}.sizes + """ + + stub: + """ + touch ${params.reference}.sizes + """ +} + +process bgcFilter { + tag "medium_mem" + publishDir "${params.outdir}/BGC", mode: "${params.publish_dir_mode}", overwrite: true + label "small" + + input: + path BED + + output: + path "${BED.simpleName}.bgc${params.bgc_threshold}.bed" + + + script: + """ + awk '\$4 >= ${params.bgc_threshold}' ${BED} > ${BED.simpleName}.bgc${params.bgc_threshold}.bed + """ + + stub: + """ + touch ${BED.simpleName}.bgc${params.bgc_threshold}.bed + """ +} + +process catsort { + 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: + """ + cat ${INFILES} | sort ${SORTEXPR} > ${OUTNAME} + """ + + stub: + """ + touch ${OUTNAME} + """ +} + +// Processes for the phastBias methods +process RENAME_BED { + + input: + path BED + path TABLE + + output: + path "${BED.simpleName}.renamed.bed" + + script: + """ + #!/usr/bin/env python + import sys + bed_fn = "${BED}" + map_fn = "${TABLE}" + + # Load the sizes and order them + conv_map = {} + with open(map_fn) as map_file: + for line in map_file: + if line[0] == "#": + continue + line = line.strip().split() + code_1 = line[0] + code_2 = line[1] + conv_map[code_1] = code_2 + + # Convert BED file + with open(bed_fn) as bed_file, open(f"${BED.simpleName}.renamed.bed", "w") as out_bed: + for line in bed_file: + line = line.strip().split() + line[0] = conv_map.get(line[0], line[0]) + out_bed.write("\\t".join(line) + "\\n") + """ +} \ No newline at end of file diff --git a/include/process/prerun.nf b/include/process/prerun.nf index 76d53cf..fb6cc7e 100755 --- a/include/process/prerun.nf +++ b/include/process/prerun.nf @@ -373,11 +373,12 @@ process make_shapeit{ process daf { tag "mutyper" label "medium" - publishDir "${params.outdir}/mutyper/daf", mode: "${params.publish_dir_mode}", overwrite: true + publishDir "${params.outdir}/${OUTDIR}/daf", mode: "${params.publish_dir_mode}", overwrite: true input: tuple val(chrom), path(vcf), path(tbi) + val OUTDIR output: path "daf.${chrom}.csv.gz" @@ -400,10 +401,11 @@ process daf { // Smile plot for the derived allele frequencies process smile { label "renv" - publishDir "${params.outdir}/mutyper/smile", mode: "${params.publish_dir_mode}", overwrite: true + publishDir "${params.outdir}/${OUTDIR}/smile", mode: "${params.publish_dir_mode}", overwrite: true input: path "dafs/*" + val OUTDIR output: path "smile.pdf" diff --git a/include/process/sdm.nf b/include/process/sdm.nf index 79c1d9f..1d9e2a2 100755 --- a/include/process/sdm.nf +++ b/include/process/sdm.nf @@ -17,7 +17,17 @@ process sdm { script: """ sdm ${vcf} ${samplelist} ${chrom} ${reffasta} sdm.${samplename}.${chrom} - gzip sdm.${samplename}.${chrom}.txt + # If a dataset is empty, we create an empty file with the appropriate name. + # This allows to check for the file size at the workflow level if the file is empty, + # and if it is, it is filtered out. If we gzip an empty file, it will still have a + # size >0 bytes, always passing the check. + if [ ! -s sdm.${samplename}.${chrom}.txt ]; then + echo "No SDM found for ${samplename} on ${chrom}" + touch sdm.${samplename}.${chrom}.txt.gz + else + gzip sdm.${samplename}.${chrom}.txt + + fi """ stub: diff --git a/include/workflow/ibd.nf b/include/workflow/ibd.nf index 829861e..8c13310 100755 --- a/include/workflow/ibd.nf +++ b/include/workflow/ibd.nf @@ -1,20 +1,13 @@ // workflow to generate the ibd windows include {get_ref_ibd; get_merge_ibd} from '../process/dependencies' -include {ibd; make_map; merge_ibd} from '../process/ibd' +include {ibd; make_map; merge_ibd; merge_ibd as merge_hbd} from '../process/ibd' +include {zcatsort as zcatsort_hbd; zcatsort as zcatsort_ibd} from '../process/ibd' workflow IBD { take: - vcf - tbi - chromosomeList + vcf_by_chr main: - // Get chromosome list - chromosomeList - .splitCsv(header: ['N','chrom']) - .map{ row-> tuple(row.N, row.chrom) } - .set{ chromosomes_ch } - // get refined-ibd if missing if (params.refinedibd){ ch_refibd = file(params.refinedibd) @@ -31,14 +24,39 @@ workflow IBD { ch_mergeibd = get_merge_ibd.out } - // make map - make_map(vcf, tbi) + // Get IBD/HBD windows + ref_ibd_ch = ibd( + make_map(vcf_by_chr), + ch_refibd + ) + + // Merge IBD windows + ibds_ch = merge_ibd( + ref_ibd_ch.ibd, + ch_mergeibd, + "IBD" + ) - // run refined-ibd - ibd(vcf, tbi, ch_refibd, chromosomes_ch) + // Merge HBD windows + hbds_ch = merge_hbd( + ref_ibd_ch.hbd, + ch_mergeibd, + "HBD" + ) - // Merge-ibd - merge_ibd(ibd.out, vcf, tbi, ch_mergeibd, make_map.out) + // Collect the HBD and IBD files + zcatsort_hbd( + hbds_ch | collect, + "HBD_regions.hbd.gz", + "IBD/", + "-k 1,1 -k2,2n" + ) + zcatsort_ibd( + ibds_ch | collect, + "IBD_regions.ibd.gz", + "IBD/", + "-k 1,1 -k2,2n" + ) emit: merge_ibd.out diff --git a/include/workflow/phylop.nf b/include/workflow/phylop.nf index 2333055..d1fa68b 100755 --- a/include/workflow/phylop.nf +++ b/include/workflow/phylop.nf @@ -1,26 +1,37 @@ include { - hal4d; - phyloFit; - phyloPtrain; + bgcFilter; + catsort as catsort_tract; + catsort as catsort_informative; + collect_maf; combine_bed; + create_maf; extract_conserved; - vcf_drop_conserved; - phyloP; - make4dmaf; - msa_view; + GENOME_INTERVALS; + hal4d; hal_genomes; - wig2bedgraph; + halSize; halTree; + make4dmaf; + msa_view as maf_to_ss; + phastBias; + phyloFit; + phyloP; + phyloPtrain; + RENAME_BED; + split_4d_bed; + vcf_drop_intervals; + wig2bedgraph; } from '../process/phylop' +include { + daf as daf_noConstrained; + daf as daf_noBGC; + smile as smile_noConstrained; + smile as smile_noBGC; +} from '../process/prerun' -workflow CONSTRAINED { - take: - vcf - tbi - chromosomeList - +workflow NEUTRAL_MODEL { main: // Load hal file if (params.hal) { @@ -38,45 +49,189 @@ workflow CONSTRAINED { } else { exit 1, "Exon BED file ${params.exon_bed} not found" } - - // Get chromosome list - chromosomes_ch = chromosomeList - .splitCsv(header: ['N','chrom']) - .map{ row-> tuple(row.N, row.chrom) } + // Hal genomes + genomes_ch = hal_genomes(hal_ch) // Extract 4d elements - maf4d = hal4d(hal_ch, exons_ch) - | combine(hal_genomes(hal_ch)) + ss4d_ch = hal4d(hal_ch, exons_ch) + | combine(genomes_ch) | map{ - hal, maf, genomes_env -> + hal, bed, genomes_env -> String genomes = genomes_env as String - [hal, maf, genomes] + [hal, genomes, bed] } + | split_4d_bed + | transpose | make4dmaf - | msa_view + | map { + _hal, genomes, maf -> + [genomes, maf] + } + | groupTuple(by: 0) + | collect_maf + | maf_to_ss // Fit model using 4D codons - model_ch = maf4d + model_ch = ss4d_ch | combine(hal_ch | halTree) | phyloFit + // + // Make maf file for BGC calculation + // + // First, create the reference sizes file + ref_sizes_ch = halSize(hal_ch) + + // then, we create the intervals and we parse them into a nested channel + intervals_ch = GENOME_INTERVALS(ref_sizes_ch) | flatten + emit: + model = model_ch + intervals = intervals_ch + sizes = ref_sizes_ch + genomes = genomes_ch +} + + +workflow CONSTRAINED { + take: + vcf_by_chr_ch + model_ch + intervals_ch + ref_sizes_ch + + main: + // Load hal file + if (params.hal) { + if (file(params.hal).exists()){ + hal_ch = Channel.fromPath(params.hal) + } else { + exit 1, "Hal file ${params.hal} not found" + } + } else { + exit 1, 'Hal file not specified!' + } + // load exon bed file + if (file(params.exon_bed).exists()){ + exons_ch = Channel.fromPath(params.exon_bed) + } else { + exit 1, "Exon BED file ${params.exon_bed} not found" + } + // Run PhyloP - phylop_ch = chromosomes_ch - | combine(hal_ch) - | combine(model_ch) - | phyloP + phylop_ch = phyloP( + hal_ch | combine(intervals_ch), + model_ch | collect + ) + | combine(ref_sizes_ch) + // Ignore empty wig files + | filter { wig, _sizes -> wig.size() > 0 } | wig2bedgraph // combine selected | collect - | combine_bed // Extract conserved - conserved_ch = phylop_ch | extract_conserved + conserved_ch = combine_bed(phylop_ch, "phylop", "PHYLOP") + | extract_conserved + + // Rename if a mapping file is provided + if (params.rename_hal_sequences){ + conversion_table_ch = file( params.rename_hal_sequences ) + conserved_ch = RENAME_BED(conserved_ch, conversion_table_ch) + } // Perform actual filtering - vcf_drop_conserved(vcf, tbi, conserved_ch) + vcf_out_ch = vcf_drop_intervals( + vcf_by_chr_ch, + conserved_ch | collect, + "non-conserved", + "PHYLOP/" + ) + + // Regenerate smile plots for the filtered VCF + smile_noConstrained( + daf_noConstrained(vcf_out_ch, "DAF/no-constrained") | collect, + "DAF/no-constrained" + ) + + emit: + vcf = vcf_out_ch +} + + +// workflow to compute the biased gene conversion (BGC) regions +workflow BGC { + take: + vcf_by_chr_ch + model_ch + intervals_ch + ref_sizes_ch + genomes_ch + + main: + // Load hal file + if (params.hal) { + if (file(params.hal).exists()){ + hal_ch = Channel.fromPath(params.hal) + } else { + exit 1, "Hal file ${params.hal} not found" + } + } else { + exit 1, 'Hal file not specified!' + } + // load exon bed file + if (file(params.exon_bed).exists()){ + exons_ch = Channel.fromPath(params.exon_bed) + } else { + exit 1, "Exon BED file ${params.exon_bed} not found" + } + maf_ch = create_maf(hal_ch | combine(intervals_ch | flatten), genomes_ch | collect) + + // Run pastBias + phastbias_ch = phastBias(maf_ch, model_ch | collect) + + // Collect bigwig outputs + large_bed_ch = phastbias_ch.wig + | combine(ref_sizes_ch) + | wig2bedgraph + | collect + // Save tract regions + tracts_bed_ch = catsort_tract ( + phastbias_ch.tracts_bed | collect, + "bgc_regions.bed", + "PHAST/gBGC/", + "-k 1,1 -k2,2n" + ) + // Save also informative regions + informative_bed_ch = catsort_informative ( + phastbias_ch.informative_bed | collect, + "bgc_informative.bed", + "PHAST/gBGC/", + "-k 1,1 -k2,2n" + ) + + // Combine all BGC regions + bgc_ch = combine_bed(large_bed_ch, "bgc", "PHASTBIAS") + + // Rename if a mapping file is provided + if (params.rename_hal_sequences){ + conversion_table_ch = file( params.rename_hal_sequences ) + tracts_bed_ch = RENAME_BED(tracts_bed_ch, conversion_table_ch) + } + + // Perform actual filtering + vcf_out_ch = vcf_drop_intervals( + vcf_by_chr_ch, + tracts_bed_ch | collect, + 'no-bgc', + "PHASTBIAS/" + ) + + // Regenerate smile plots for the filtered VCF + smile_noBGC( + daf_noBGC(vcf_out_ch, "DAF/no-BGC") | collect, + "DAF/no-BGC" + ) emit: - vcf = vcf_drop_conserved.out.vcf - tbi = vcf_drop_conserved.out.tbi + vcf = vcf_out_ch } \ No newline at end of file diff --git a/include/workflow/sdm.nf b/include/workflow/sdm.nf index a5aae7b..ae0a26d 100755 --- a/include/workflow/sdm.nf +++ b/include/workflow/sdm.nf @@ -38,6 +38,7 @@ workflow SDM { raw_sdm = breeds_ch | combine( sdm( combined_ch, reffasta, reffai ) + | filter { _group, sdm -> sdm.size() > 0 } // Drop empty files; if all files for a breed are empty, the breed will be dropped from the pipeline | groupTuple(by: 0), by: 0 ) diff --git a/include/workflow/vcfpreprocess.nf b/include/workflow/vcfpreprocess.nf index e949243..420968d 100755 --- a/include/workflow/vcfpreprocess.nf +++ b/include/workflow/vcfpreprocess.nf @@ -95,8 +95,11 @@ workflow PREPROCESS { ) // Generate derived allele frequency plots prior DAF filtering. - ancestral_ch | daf | collect | smile - + smile( + daf(ancestral_ch, "DAF/filtered") | collect, + "DAF/filtered" + ) + // Filter by derived allele freq. filtered_ch = ancestral_ch | daf_filter diff --git a/main.nf b/main.nf index 5cdba53..25a226c 100755 --- a/main.nf +++ b/main.nf @@ -12,7 +12,7 @@ include { RELATE } from './include/workflow/relate' include { SDM } from './include/workflow/sdm' include { GONE } from './include/workflow/gone' include { IBD } from './include/workflow/ibd' -include { CONSTRAINED } from './include/workflow/phylop' +include { BGC; CONSTRAINED; NEUTRAL_MODEL } from './include/workflow/phylop' include { get_masks } from './include/process/prerun' include { chromosomeList } from './include/process/prerun' @@ -42,31 +42,31 @@ workflow { log.info """\ Nextflow Mutation Spectra v${workflow.manifest.version} ================================================================================= - variants : $params.variants - idx : $params.idx - output folder : $params.outdir - hal : $params.hal - reference : $params.reference - target : $params.target - mutyper : $params.mutyper - sdm : $params.sdm - relate : $params.relate - species : $params.species - k : $params.k - Ne subset : $params.ne_subset - Intergen. time : $params.intergen_time - Mut. rate : $params.mutation_rate - Min. pop. size : $params.min_pop_size - imputation sfw : $params.imputation - coding : $params.coding - noncoding : $params.noncoding - annotation : $params.annotation - pops_folder : $params.pops_folder - pop labels file : $params.poplabels - chromosome list : $params.chr_list - cactus url : $params.cactus_url - exons : $params.exon_bed - constrained : $params.constrained + variants : $params.variants + idx : $params.idx + output folder : $params.outdir + hal : $params.hal + reference : $params.reference + target : $params.target + mutyper : $params.mutyper + sdm : $params.sdm + relate : $params.relate + species : $params.species + k : $params.k + Ne subset : $params.ne_subset + Intergen. time : $params.intergen_time + Mut. rate : $params.mutation_rate + Min. pop. size : $params.min_pop_size + imputation sfw : $params.imputation + coding : $params.coding + noncoding : $params.noncoding + annotation : $params.annotation + pops_folder : $params.pops_folder + pop labels file : $params.poplabels + chromosome list : $params.chr_list + cactus url : $params.cactus_url + exons : $params.exon_bed + remove_constrained : $params.remove_constrained relate dir : $params.relate""" if (params.neval){ log.info """Ne value : $params.neval""" @@ -122,7 +122,8 @@ workflow { } if (!params.hal) { exit 1, 'Hal file not specified and ancestral not specified!' } if (!params.hal && !params.reference_fna && !params.ancestral_fna) { exit 1, 'Ancestral and reference genomes not specified!' } - if (params.constrained && !params.exon_bed) { exit 1, 'Requested hal4d algorithm, but no bed with exons specified!' } + if (params.remove_constrained && !params.exon_bed) { exit 1, 'Requested constrained elements removal sub-workflow, but no bed with exons specified!' } + if (params.bgc && !params.exon_bed) { exit 1, 'Requested BGC sub-workflow, but no bed with exons specified!' } if (!params.species){ throw new Exception("Parameter --species is required for file naming.") } // Generate the ancestral fasta @@ -148,30 +149,46 @@ workflow { ch_var_idx_new = PREPROCESS.out.tbi ch_chr_lists = PREPROCESS.out.chroms vcf_by_chr = PREPROCESS.out.vcf_by_chr - vcf_chunks_ch = PREPROCESS.out.chunks_ch - - // Get constrined elements and remove variants in them - if ( params.constrained ){ - CONSTRAINED(ch_var_new, ch_var_idx_new, ch_chr_lists) - ch_var_new = CONSTRAINED.out[0] - ch_var_idx_new = CONSTRAINED.out[1] - } - - // Generate IBDs if requested - if (params.compute_ibd){ - IBD(ch_var_new, ch_var_idx_new, ch_chr_lists) - } } else { ch_var_new = ch_var - ch_var_idx_new = ch_var_idx - chromosomeList( ch_var, ch_var_idx ) - ch_chr_lists = chromosomeList.out + ch_var_idx_new = ch_var_idx + ch_chr_lists = chromosomeList( ch_var, ch_var_idx ) + chromosomes_ch = ch_chr_lists + | splitCsv(header: ['N','chrom']) + | map{ row-> [row.chrom] } + vcf_by_chr = chromosomes_ch | combine( ch_var_new | combine( ch_var_idx ) ) } + // Define neutral model once for both constrained and gBGC workflows + if ( params.remove_constrained || params.bgc ){ + neutral_model = NEUTRAL_MODEL() + model_ch = neutral_model.model + intervals_ch = neutral_model.intervals + ref_sizes_ch = neutral_model.sizes + genomes_ch = neutral_model.genomes + } + + // Get constrined elements and remove variants in them + if ( params.remove_constrained ){ + phastcons_ch = CONSTRAINED(vcf_by_chr, model_ch, intervals_ch, ref_sizes_ch) + vcf_by_chr = phastcons_ch.vcf + } + + // Get variants not in BGC candidate regions + if ( params.bgc ){ + bgc_ch = BGC(vcf_by_chr, model_ch, intervals_ch, ref_sizes_ch, genomes_ch) + vcf_by_chr = bgc_ch.vcf + } + + // Generate IBDs if requested + if (params.compute_ibd){ + IBD(vcf_by_chr) + } + // Run GONE to calculate Ne, if requested if (params.gone){ GONE( ch_var_new, ch_var_idx_new ) diff --git a/nextflow.config b/nextflow.config index eb58140..da8513c 100755 --- a/nextflow.config +++ b/nextflow.config @@ -23,10 +23,13 @@ params { ancestral_fna = false ancestral_only = false greedy = true + cactus_hal2maf = false + // Spectrum mutyper = true sdm = false relate = false + // Preprocessing vep = true vcf_is_filtered = false @@ -94,13 +97,20 @@ params { vep_args = null // Beagle settings - refined_ibd_params = 'window=10 trim=0.015 length=0.2 lod=4' - merge_ibd_params = '1 0.6' + recombination_rate = 1e-8 + refined_ibd_window = 0.5 + refined_ibd_trim = 0.005 + refined_ibd_length = 0.1 + refined_ibd_lod = 4 + merge_ibd_gap = 0.25 + merge_ibd_discordance = 10 cactus_version = null // Conserved elements detection - constrained = false - exon_bed = null + remove_constrained = false + bgc = false + bgc_threshold = 0.5 + exon_bed = null // Other mamba = false @@ -127,7 +137,7 @@ manifest { description = 'Nextflow mutation spectra analysis workflow.' mainScript = 'main.nf' nextflowVersion = '>=21.04.0' - version = '1.1.6a' + version = '1.2.0' } profiles { diff --git a/singularity.def b/singularity.def index a074355..3f7c795 100644 --- a/singularity.def +++ b/singularity.def @@ -18,13 +18,12 @@ mamba install -n nspectra -y \ bioconda::bedtools bioconda::plink=1.90 mamba install -n nspectra -y bioconda::ucsc-twobitinfo bioconda::ucsc-fatotwobit \ bioconda::ucsc-wigtobigwig bioconda::ucsc-bigwigtobedgraph -mamba install -n nspectra -y bioconda::vcflib==1.0.3 mamba install -n nspectra -y bioconda::shapeit4 mamba install -n nspectra -y bioconda::shapeit5 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 mamba install -n nspectra -y \ conda-forge::r-base>=4.1.0 \ conda-forge::r-cowplot \ @@ -48,32 +47,40 @@ rm /tmp/env.tar # so now fix up paths: /venv/bin/conda-unpack +# The cargo compiler image; we can use alpine as it is small. +Bootstrap: docker +From: rust:1.67-alpine +Stage: cargo_build +WORKDIR /opt/ + +%post +cd /opt +apk add --no-cache git +git clone https://github.com/julibeg/smakcr +cd smakcr && \ + cargo build --release && \ + cp target/release/smakcr /opt # The runtime-stage image; we can use Debian as the # base image since the Conda env also includes Python # for us. - - Bootstrap: docker From: ubuntu:24.04 Stage: runtime -%files from build -/venv /venv %post - # Install procps in debian to make it compatible with reporting apt-get update && apt install -y file jellyfish procps g++ curl git wget parallel && apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* curl -o /usr/local/bin/datasets https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/LATEST/linux-amd64/datasets && \ chmod a+x /usr/local/bin/datasets -# Copy /venv from the previous stage: +%files from build +/venv /venv + +%files from cargo_build +/opt/smakcr /usr/local/bin/smakcr -# When image is run, run the code with the environment -# activated: -PATH=/venv/bin/:$PATH -SHELL ["/bin/bash", "-c"] %environment export PATH=/venv/bin/:$PATH %runscript diff --git a/test/demo.hal b/test/demo.hal index ddcb9d6..d2722f6 100755 Binary files a/test/demo.hal and b/test/demo.hal differ