Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-27T22:23:26.280Z
Answers.LastActivityDate: 2018-03-27T22:23:26.280Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-21T15:22:24.297Z
Answers.LastActivityDate: 2018-03-21T16:24:05.467Z
Answers.LastEditDate: 2018-03-21T16:24:05.467Z
Answers.UserDisplayName: bli
Body: <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>
  
Tags: <bed><sam><text-processing>
Title: How can I easily get the read size distribution of reads mapping on a certain set of regions?
ThreadText: # 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](https://www.biostars.org/p/98407/#98408) 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](https://github.com/bedops/bedops) `bedmap` and `bam2bed`, along with the `hist` functionality in [`bashplotlib`](https://github.com/glamp/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    |
  ------------------------------
  
  ```
  
  
  ---
  
  