<p>Via <a href="https://github.com/bedops/bedops" rel="nofollow noreferrer">BEDOPS</a> <code>bedmap</code> and <code>bam2bed</code>, along with the <code>hist</code> functionality in <a href="https://github.com/glamp/bashplotlib" rel="nofollow noreferrer"><code>bashplotlib</code></a></p>

<pre><code>$ bedmap --echo-map-size --multidelim '\n' intervals.bed &lt;(bam2bed &lt; reads.bam) | hist

 77|           o
 73|           o
 69|           o
 65|           o
 61|           o
 57|           o
 53|           o
 49|           o
 45|           o
 41|           o
 37|           o
 33|           o
 29|           o
 25|           o
 21|           o
 17|          oo
 13|         ooo
  9|         ooo
  5|         ooo
  1| o    o oooo
    -----------

------------------------------
|          Summary           |
------------------------------
|     observations: 118      |
|    min value: 16.000000    |
|      mean : 34.415254      |
|    max value: 36.000000    |
------------------------------
</code></pre>

CC BY-SA 3.0
2018-03-27T22:23:26.280Z
2018-03-27T22:23:26.280Z
Alex Reynolds
<p>This can be done rather simply by combining <a href="https://www.biostars.org/p/98407/#98408" rel="nofollow noreferrer">this samtools based answer</a> and bioawk:</p>

<pre><code>samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1
</code></pre>

<p>The <code>samtools</code> command extracts the alignements in the region given with the <code>-L</code> option, in sam format. The <code>-c sam</code> option of <code>bioawk</code> makes the <code>$seq</code> variable (corresponding to the read sequence field) available in the awk command.</p>

<p>The final <code>sort</code> ensures the size classes are ordered numerically.</p>

<p>We obtain a two-column <code>size    count</code> table.</p>

<p>Examples with timings:</p>

<pre><code>$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s
</code></pre>

<p>Since sam format is actually simple to parse, one can try to use <code>mawk</code> instead for slightly higher speed (the sequence being in the 10th field):</p>

<pre><code>samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1
</code></pre>

<p>(Contrary to bioawk, the default output field separator is not a tabulation.)</p>

<p>Example:</p>

<pre><code>$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s
</code></pre>

CC BY-SA 3.0
2018-03-21T15:22:24.297Z
2018-03-21T16:24:05.467Z
2018-03-21T16:24:05.467Z
bli
<p>Suppose I have a BAM file indicating where reads in a library have mapped, and a bed file describing a set of genomic regions.</p>

<p>Is there a way to easily get the size distribution of the reads mapping on this set of regions?</p>

<bed><sam><text-processing>
How can I easily get the read size distribution of reads mapping on a certain set of regions?
# How can I easily get the read size distribution of reads mapping on a certain set of regions?

(Asked by: username_0 on 2018-03-21)

Suppose I have a BAM file indicating where reads in a library have mapped, and a bed file describing a set of genomic regions.

Is there a way to easily get the size distribution of the reads mapping on this set of regions?

**Tags:** bed, sam, text-processing

---

## 2 answers

### Answer by username_0 on 2018-03-21. Score: 4 | Accepted answer

This can be done rather simply by combining this samtools based answer and bioawk:
samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1

The samtools command extracts the alignements in the region given with the -L option, in sam format. The -c sam option of bioawk makes the $seq variable (corresponding to the read sequence field) available in the awk command.
The final sort ensures the size classes are ordered numerically.
We obtain a two-column size    count table.
Examples with timings:
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s

Since sam format is actually simple to parse, one can try to use mawk instead for slightly higher speed (the sequence being in the 10th field):
samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1

(Contrary to bioawk, the default output field separator is not a tabulation.)
Example:
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s


---

### Answer by username_1 on 2018-03-27. Score: 1

Via BEDOPS bedmap and bam2bed, along with the hist functionality in bashplotlib
$ bedmap --echo-map-size --multidelim '\n' intervals.bed <(bam2bed < reads.bam) | hist

 77|           o
 73|           o
 69|           o
 65|           o
 61|           o
 57|           o
 53|           o
 49|           o
 45|           o
 41|           o
 37|           o
 33|           o
 29|           o
 25|           o
 21|           o
 17|          oo
 13|         ooo
  9|         ooo
  5|         ooo
  1| o    o oooo
    -----------

------------------------------
|          Summary           |
------------------------------
|     observations: 118      |
|    min value: 16.000000    |
|      mean : 34.415254      |
|    max value: 36.000000    |
------------------------------


---

