diff --git a/CHANGELOG.md b/CHANGELOG.md index 3e9f40fc..8c1af805 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,6 +29,7 @@ * `bedtools`: - `bedtools/bedtools_intersect`: Allows one to screen for overlaps between two sets of genomic features (PR #94). - `bedtools/bedtools_sort`: Sorts a feature file (bed/gff/vcf) by chromosome and other criteria (PR #98). + - `bedtools/bedtools_merge`: Merges overlapping BED/GFF/VCF entries into a single interval (PR #118). - `bedtools/bedtools_bamtofastq`: Convert BAM alignments to FASTQ files (PR #101). - `bedtools/bedtools_bedtobam`: Converts genomic feature records (bed/gff/vcf) to BAM format (PR #111). diff --git a/src/bedtools/bedtools_merge/config.vsh.yaml b/src/bedtools/bedtools_merge/config.vsh.yaml new file mode 100644 index 00000000..45e4a01d --- /dev/null +++ b/src/bedtools/bedtools_merge/config.vsh.yaml @@ -0,0 +1,160 @@ +name: bedtools_merge +namespace: bedtools +description: | + Merges overlapping BED/GFF/VCF entries into a single interval. +links: + documentation: https://bedtools.readthedocs.io/en/latest/content/tools/merge.html + repository: https://github.com/arq5x/bedtools2 + homepage: https://bedtools.readthedocs.io/en/latest/# + issue_tracker: https://github.com/arq5x/bedtools2/issues +references: + doi: 10.1093/bioinformatics/btq033 +license: MIT +requirements: + commands: [bedtools] +authors: + - __merge__: /src/_authors/theodoro_gasperin.yaml + roles: [ author, maintainer ] + +argument_groups: + - name: Inputs + arguments: + - name: --input + alternatives: -i + type: file + description: Input file (BED/GFF/VCF) to be merged. + required: true + + - name: Outputs + arguments: + - name: --output + type: file + direction: output + description: Output merged file BED to be written. + required: true + + - name: Options + arguments: + - name: --strand + alternatives: -s + type: boolean_true + description: | + Force strandedness. That is, only merge features + that are on the same strand. + - By default, merging is done without respect to strand. + + - name: --specific_strand + alternatives: -S + type: string + choices: ["+", "-"] + description: | + Force merge for one specific strand only. + Follow with + or - to force merge from only + the forward or reverse strand, respectively. + - By default, merging is done without respect to strand. + + - name: --distance + alternatives: -d + type: integer + description: | + Maximum distance between features allowed for features + to be merged. + - Def. 0. That is, overlapping & book-ended features are merged. + - (INTEGER) + - Note: negative values enforce the number of b.p. required for overlap. + + - name: --columns + alternatives: -c + type: integer + description: | + Specify columns from the B file to map onto intervals in A. + Default: 5. + Multiple columns can be specified in a comma-delimited list. + + - name: --operation + alternatives: -o + type: string + description: | + Specify the operation that should be applied to -c. + Valid operations: + sum, min, max, absmin, absmax, + mean, median, mode, antimode + stdev, sstdev + collapse (i.e., print a delimited list (duplicates allowed)), + distinct (i.e., print a delimited list (NO duplicates allowed)), + distinct_sort_num (as distinct, sorted numerically, ascending), + distinct_sort_num_desc (as distinct, sorted numerically, desscending), + distinct_only (delimited list of only unique values), + count + count_distinct (i.e., a count of the unique values in the column), + first (i.e., just the first value in the column), + last (i.e., just the last value in the column), + Default: sum + Multiple operations can be specified in a comma-delimited list. + + If there is only column, but multiple operations, all operations will be + applied on that column. Likewise, if there is only one operation, but + multiple columns, that operation will be applied to all columns. + Otherwise, the number of columns must match the the number of operations, + and will be applied in respective order. + E.g., "-c 5,4,6 -o sum,mean,count" will give the sum of column 5, + the mean of column 4, and the count of column 6. + The order of output columns will match the ordering given in the command. + + - name: --delimiter + alternatives: -delim + type: string + description: | + Specify a custom delimiter for the collapse operations. + example: "|" + default: "," + + - name: --precision + alternatives: -prec + type: integer + description: | + Sets the decimal precision for output (Default: 5). + + - name: --bed + type: boolean_true + description: | + If using BAM input, write output as BED. + + - name: --header + type: boolean_true + description: | + Print the header from the A file prior to results. + + - name: --no_buffer + alternatives: -nobuf + type: boolean_true + description: | + Disable buffered output. Using this option will cause each line + of output to be printed as it is generated, rather than saved + in a buffer. This will make printing large output files + noticeably slower, but can be useful in conjunction with + other software tools and scripts that need to process one + line of bedtools output at a time. + +resources: + - type: bash_script + path: script.sh + +test_resources: + - type: bash_script + path: test.sh + - path: test_data + +engines: + - type: docker + image: debian:stable-slim + setup: + - type: apt + packages: [bedtools, procps] + - type: docker + run: | + echo "bedtools: \"$(bedtools --version | sed -n 's/^bedtools //p')\"" > /var/software_versions.txt + +runners: + - type: executable + - type: nextflow \ No newline at end of file diff --git a/src/bedtools/bedtools_merge/help.txt b/src/bedtools/bedtools_merge/help.txt new file mode 100644 index 00000000..bc78fc67 --- /dev/null +++ b/src/bedtools/bedtools_merge/help.txt @@ -0,0 +1,85 @@ +```bash +bedtools merge +``` + +Tool: bedtools merge (aka mergeBed) +Version: v2.30.0 +Summary: Merges overlapping BED/GFF/VCF entries into a single interval. + +Usage: bedtools merge [OPTIONS] -i + +Options: + -s Force strandedness. That is, only merge features + that are on the same strand. + - By default, merging is done without respect to strand. + + -S Force merge for one specific strand only. + Follow with + or - to force merge from only + the forward or reverse strand, respectively. + - By default, merging is done without respect to strand. + + -d Maximum distance between features allowed for features + to be merged. + - Def. 0. That is, overlapping & book-ended features are merged. + - (INTEGER) + - Note: negative values enforce the number of b.p. required for overlap. + + -c Specify columns from the B file to map onto intervals in A. + Default: 5. + Multiple columns can be specified in a comma-delimited list. + + -o Specify the operation that should be applied to -c. + Valid operations: + sum, min, max, absmin, absmax, + mean, median, mode, antimode + stdev, sstdev + collapse (i.e., print a delimited list (duplicates allowed)), + distinct (i.e., print a delimited list (NO duplicates allowed)), + distinct_sort_num (as distinct, sorted numerically, ascending), + distinct_sort_num_desc (as distinct, sorted numerically, desscending), + distinct_only (delimited list of only unique values), + count + count_distinct (i.e., a count of the unique values in the column), + first (i.e., just the first value in the column), + last (i.e., just the last value in the column), + Default: sum + Multiple operations can be specified in a comma-delimited list. + + If there is only column, but multiple operations, all operations will be + applied on that column. Likewise, if there is only one operation, but + multiple columns, that operation will be applied to all columns. + Otherwise, the number of columns must match the the number of operations, + and will be applied in respective order. + E.g., "-c 5,4,6 -o sum,mean,count" will give the sum of column 5, + the mean of column 4, and the count of column 6. + The order of output columns will match the ordering given in the command. + + + -delim Specify a custom delimiter for the collapse operations. + - Example: -delim "|" + - Default: ",". + + -prec Sets the decimal precision for output (Default: 5) + + -bed If using BAM input, write output as BED. + + -header Print the header from the A file prior to results. + + -nobuf Disable buffered output. Using this option will cause each line + of output to be printed as it is generated, rather than saved + in a buffer. This will make printing large output files + noticeably slower, but can be useful in conjunction with + other software tools and scripts that need to process one + line of bedtools output at a time. + + -iobuf Specify amount of memory to use for input buffer. + Takes an integer argument. Optional suffixes K/M/G supported. + Note: currently has no effect with compressed files. + +Notes: + (1) The input file (-i) file must be sorted by chrom, then start. + + + + +***** ERROR: No input file given. Exiting. ***** diff --git a/src/bedtools/bedtools_merge/script.sh b/src/bedtools/bedtools_merge/script.sh new file mode 100644 index 00000000..db50dd83 --- /dev/null +++ b/src/bedtools/bedtools_merge/script.sh @@ -0,0 +1,35 @@ +#!/bin/bash + +## VIASH START +## VIASH END + +# Exit on error +set -eo pipefail + +# Unset parameters +unset_if_false=( + par_strand + par_bed + par_header + par_no_buffer +) + +for par in ${unset_if_false[@]}; do + test_val="${!par}" + [[ "$test_val" == "false" ]] && unset $par +done + +# Execute bedtools merge with the provided arguments +bedtools merge \ + ${par_strand:+-s} \ + ${par_specific_strand:+-S "$par_specific_strand"} \ + ${par_bed:+-bed} \ + ${par_header:+-header} \ + ${par_no_buffer:+-nobuf} \ + ${par_distance:+-d "$par_distance"} \ + ${par_columns:+-c "$par_columns"} \ + ${par_operation:+-o "$par_operation"} \ + ${par_delimiter:+-delim "$par_delimiter"} \ + ${par_precision:+-prec "$par_precision"} \ + -i "$par_input" \ + > "$par_output" diff --git a/src/bedtools/bedtools_merge/test.sh b/src/bedtools/bedtools_merge/test.sh new file mode 100644 index 00000000..e2b46c15 --- /dev/null +++ b/src/bedtools/bedtools_merge/test.sh @@ -0,0 +1,222 @@ +#!/bin/bash + +# exit on error +set -eo pipefail + +## VIASH START +meta_executable="target/executable/bedtools/bedtools_sort/bedtools_merge" +meta_resources_dir="src/bedtools/bedtools_merge" +## VIASH END + +# directory of the bam file +test_data="$meta_resources_dir/test_data" + +############################################# +# helper functions +assert_file_exists() { + [ -f "$1" ] || { echo "File '$1' does not exist" && exit 1; } +} +assert_file_not_empty() { + [ -s "$1" ] || { echo "File '$1' is empty but shouldn't be" && exit 1; } +} +assert_file_contains() { + grep -q "$2" "$1" || { echo "File '$1' does not contain '$2'" && exit 1; } +} +assert_identical_content() { + diff -a "$2" "$1" \ + || (echo "Files are not identical!" && exit 1) +} +############################################# + +# Create directories for tests +echo "Creating Test Data..." +TMPDIR=$(mktemp -d "$meta_temp_dir/XXXXXX") +function clean_up { + [[ -d "$TMPDIR" ]] && rm -r "$TMPDIR" +} +trap clean_up EXIT + +# Create and populate example files +printf "chr1\t100\t200\nchr1\t150\t250\nchr1\t300\t400\n" > "$TMPDIR/featureA.bed" +printf "chr1\t100\t200\ta1\t1\t+\nchr1\t180\t250\ta2\t2\t+\nchr1\t250\t500\ta3\t3\t-\nchr1\t501\t1000\ta4\t4\t+\n" > "$TMPDIR/featureB.bed" +printf "chr1\t100\t200\ta1\t1.9\t+\nchr1\t180\t250\ta2\t2.5\t+\nchr1\t250\t500\ta3\t3.3\t-\nchr1\t501\t1000\ta4\t4\t+\n" > "$TMPDIR/feature_precision.bed" + +# Create and populate feature.gff file +printf "##gff-version 3\n" > "$TMPDIR/feature.gff" +printf "chr1\t.\tgene\t1000\t2000\t.\t+\t.\tID=gene1;Name=Gene1\n" >> "$TMPDIR/feature.gff" +printf "chr1\t.\texon\t1000\t1200\t.\t+\t.\tID=exon1;Parent=transcript1\n" >> "$TMPDIR/feature.gff" +printf "chr1\t.\tCDS\t1000\t1200\t.\t+\t0\tID=cds1;Parent=transcript1\n" >> "$TMPDIR/feature.gff" +printf "chr1\t.\tCDS\t1500\t1700\t.\t+\t2\tID=cds2;Parent=transcript1\n" >> "$TMPDIR/feature.gff" +printf "chr2\t.\texon\t1500\t1700\t.\t+\t.\tID=exon2;Parent=transcript1\n" >> "$TMPDIR/feature.gff" +printf "chr3\t.\tmRNA\t1000\t2000\t.\t+\t.\tID=transcript1;Parent=gene1\n" >> "$TMPDIR/feature.gff" + +# Create expected output files +printf "chr1\t100\t250\nchr1\t300\t400\n" > "$TMPDIR/expected.bed" +printf "chr1\t100\t250\nchr1\t250\t500\nchr1\t501\t1000\n" > "$TMPDIR/expected_strand.bed" +printf "chr1\t100\t250\nchr1\t501\t1000\n" > "$TMPDIR/expected_specific_strand.bed" +printf "chr1\t128\t228\nchr1\t428\t528\n" > "$TMPDIR/expected_bam.bed" +printf "chr1\t100\t400\n" > "$TMPDIR/expected_distance.bed" +printf "chr1\t100\t500\t2\t1\t3\nchr1\t501\t1000\t4\t4\t4\n" > "$TMPDIR/expected_operation.bed" +printf "chr1\t100\t500\ta1|a2|a3\nchr1\t501\t1000\ta4\n" > "$TMPDIR/expected_delim.bed" +printf "chr1\t100\t500\t2.567\nchr1\t501\t1000\t4\n" > "$TMPDIR/expected_precision.bed" +printf "##gff-version 3\nchr1\t999\t2000\nchr2\t1499\t1700\nchr3\t999\t2000\n" > "$TMPDIR/expected_header.bed" + +# Test 1: Default sort on BED file +mkdir "$TMPDIR/test1" && pushd "$TMPDIR/test1" > /dev/null + +echo "> Run bedtools_merge on BED file" +"$meta_executable" \ + --input "../featureA.bed" \ + --output "output.bed" + +# # checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected.bed" +echo "- test1 succeeded -" + +popd > /dev/null + +# Test 2: strand option +mkdir "$TMPDIR/test2" && pushd "$TMPDIR/test2" > /dev/null + +echo "> Run bedtools_merge on BED file with strand option" +"$meta_executable" \ + --input "../featureB.bed" \ + --output "output.bed" \ + --strand + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected_strand.bed" +echo "- test2 succeeded -" + +popd > /dev/null + +# Test 3: specific strand option +mkdir "$TMPDIR/test3" && pushd "$TMPDIR/test3" > /dev/null + +echo "> Run bedtools_merge on BED file with specific strand option" +"$meta_executable" \ + --input "../featureB.bed" \ + --output "output.bed" \ + --specific_strand "+" + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected_specific_strand.bed" +echo "- test3 succeeded -" + +popd > /dev/null + +# Test 4: BED option +mkdir "$TMPDIR/test4" && pushd "$TMPDIR/test4" > /dev/null + +echo "> Run bedtools_merge on BAM file with BED option" +"$meta_executable" \ + --input "$test_data/feature.bam" \ + --output "output.bed" \ + --bed + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected_bam.bed" +echo "- test4 succeeded -" + +popd > /dev/null + +# Test 5: distance option +mkdir "$TMPDIR/test5" && pushd "$TMPDIR/test5" > /dev/null + +echo "> Run bedtools_merge on BED file with distance option" +"$meta_executable" \ + --input "../featureA.bed" \ + --output "output.bed" \ + --distance -5 + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected.bed" +echo "- test5 succeeded -" + +popd > /dev/null + +# Test 6: columns option & operation option +mkdir "$TMPDIR/test6" && pushd "$TMPDIR/test6" > /dev/null + +echo "> Run bedtools_merge on BED file with columns & operation options" +"$meta_executable" \ + --input "../featureB.bed" \ + --output "output.bed" \ + --columns 5 \ + --operation "mean,min,max" + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected_operation.bed" +echo "- test6 succeeded -" + +popd > /dev/null + +# Test 7: delimeter option +mkdir "$TMPDIR/test7" && pushd "$TMPDIR/test7" > /dev/null + +echo "> Run bedtools_merge on BED file with delimeter option" +"$meta_executable" \ + --input "../featureB.bed" \ + --output "output.bed" \ + --columns 4 \ + --operation "collapse" \ + --delimiter "|" + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected_delim.bed" +echo "- test7 succeeded -" + +popd > /dev/null + +# Test 8: precision option +mkdir "$TMPDIR/test8" && pushd "$TMPDIR/test8" > /dev/null + +echo "> Run bedtools_merge on BED file with precision option" +"$meta_executable" \ + --input "../feature_precision.bed" \ + --output "output.bed" \ + --columns 5 \ + --operation "mean" \ + --precision 4 + +# checks +assert_file_exists "output.bed" +assert_file_not_empty "output.bed" +assert_identical_content "output.bed" "../expected_precision.bed" +echo "- test8 succeeded -" + +popd > /dev/null + +# Test 9: header option +mkdir "$TMPDIR/test9" && pushd "$TMPDIR/test9" > /dev/null + +echo "> Run bedtools_merge on GFF file with header option" +"$meta_executable" \ + --input "../feature.gff" \ + --output "output.gff" \ + --header + +# checks +assert_file_exists "output.gff" +assert_file_not_empty "output.gff" +assert_identical_content "output.gff" "../expected_header.bed" +echo "- test9 succeeded -" + +popd > /dev/null + +echo "---- All tests succeeded! ----" +exit 0 diff --git a/src/bedtools/bedtools_merge/test_data/feature.bam b/src/bedtools/bedtools_merge/test_data/feature.bam new file mode 100644 index 00000000..3d56a631 Binary files /dev/null and b/src/bedtools/bedtools_merge/test_data/feature.bam differ