This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Bedtools: merge all features with the same id

I would like to use bedtools merge to collapse together all the features sharing a same gene_id in my bed file (which contains the annotation of various genes - the name column (4th one) also corresponds to the gene_id). Due to splicing, elements may be quite distant.... ​

chr14    49894259    49895806    ENSMUST00000053290    0.000000    -    mm10_ensGene    exon    .    gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14    49894873    49894876    ENSMUST00000053290    0.000000    -    mm10_ensGene    stop_codon    .    gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14    49894876    49895800    ENSMUST00000053290    0.000000    -    mm10_ensGene    CDS    0    gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14    49895797    49895800    ENSMUST00000053290    0.000000    -    mm10_ensGene    start_codon    .    gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14    49901908    49901941    ENSMUST00000053290    0.000000    -    mm10_ensGene    exon    .    gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
annotation bedtools bed
• 13,728 views • link
updated by asalimih • asalimih • written by cyril-cros cyril-cros

Are you starting out with a GTF or a bed file? I think I have done something similar starting out with a GTF and using bedtools merge.

• 0 views • link
written by A. Domingues A. Domingues

What are you really trying to achieve in the end? For example, do you need coding gene lengths because you want to calculate transcript-size corrected gene expression values from RNA sequencing (aka [FRC]PKM?) In that case, have a look at this post.

• 0 views • link
updated by Ram Ram • written by Irsan Irsan

Nice link. I wanted to collapse my annotation to a more convenient form to intersect a de novo annotation to a reference annotation with gene_ids (see How To Get Annotation For Bed File From Another Bed File). I do this to assign the official Ensembl ids to my new annotation. I don't really need the whole transcript...

However, some parts of my annotation have a rather insufficient coverage - your script will help me a lot to identify them.

• 0 views • link
updated by Ram Ram • written by cyril-cros cyril-cros

2 answers

Hi- You don't say what to do with the collapsed features. This piece of code collapses the gene ids and counts how many features have been grouped together along with the min and max coordinates. Note it uses bedtools groupby not merge.

Sample data:

cat genes.bed
chr14    49894259    49895806    ENSMUST00000053290    0.000000    ...
chr14    49894873    49894876    ENSMUST00000053290    0.000000    ...
chr14    49894876    49895800    ENSMUST00000053291    0.000000    ...
chr14    49895797    49895800    ENSMUST00000053291    0.000000    ...
chr14    49901908    49901941    ENSMUST00000053291    0.000000    ...

Example output:

sort -k4,4 genes.bed \
| groupBy -g 1,4 -c 4,2,3 -o count,min,max \
| awk -v OFS='\t' '{print $1, $4, $5, $2, $3}'

chr14    49894259    49895806    ENSMUST00000053290    2
chr14    49894876    49901941    ENSMUST00000053291    3
• 0 views • link
updated by Ram Ram • written by dariober dariober

Thanks, it does the job nicely.

• 0 views • link
written by cyril-cros cyril-cros

Thanks! That is exactly what I was looking for. Best

• 1 views • link
written by backfish • backfish

The code didn't work for me at first. The groupBy function summarized all lines of the bed file and the output file contained one entry which covered the whole genome. It took me a while to figure out that the problem was the recent bedtools version 2.26.0. When I switched back to version 2.25.0 everything worked as expected. Thanks for your solution anyway, it helped me a lot.

• 0 views • link
written by akirbis • akirbis

How is it possible to run this code on newer versions of bedtools without summarizing all lines of each ID to one?

• 0 views • link
written by asalimih • asalimih

Here's an untested way to do this with BEDOPS. Merging intervals in a BED file is "destructive" to input regions, so there's no need to parse out anything beyond the first four fields.

#!/bin/bash +x

inFn=$1
outDir=$2
tmpDir="/tmp/split"
mkdir -p $tmpDir
mkdir -p $outDir
awk -v tmpDir=$tmpDir '
    BEGIN {
        while (getline < "'"$inFn"'") {
            n = split($0, elems, "\t");
            chr = elems[1];
            start = elems[2];
            stop = elems[3];
            id = elems[4];
            out_fn = tmpDir"/"id".id.bed";
            print chr"\t"start"\t"stop >> out_fn;
        }
        close("'"$inFn"'");
    }
'
for $perIdFn in `ls $tmpDir/*.id.bed`;
do
    perIdBaseFn=$(basename "$perIdFn")
    bedops --merge $perIdFn > $outDir/$perIdBaseFn.merged.bed
done
rm -rf $tmpDir

Usage:

$ ./splitByIdAndMerge.bash unmerged_features.bed merged_features_directory
• 0 views • link
updated by Ram Ram • written by Alex Reynolds Alex Reynolds

Thanks for your help. My main issue was that I could not easily obtain a relation in the style of one gene, one large interval. I wanted to find which genes were included in a large bed file.

By only taking the start_codon feature, I have a unique location for each gene. Granted, I miss any small overlap... I will still try your solution and validate it later if it works. A simple python parser might also work if I sort by gene_id my entries...

• 0 views • link
updated by Ram Ram • written by cyril-cros cyril-cros

You might look into BEDOPS bedmap. Given an interval and a set of features (both sorted per BEDOPS sort-bed), bedmap will report all the overlapping features for a given interval. You can then chop up results with Python, Perl, etc.

• 0 views • link
updated by Ram Ram • written by Alex Reynolds Alex Reynolds

Log in to answer this question.

More posts like this