Answers.Body: <p>Your file's headers are also being passed through <code>awk</code>. Header lines' fields are all fairly short (and these lines may have fewer than ten fields, so <code>length($10)</code> would be zero), so they are all being removed by the awk conditions of the form <code>length($10) &gt;= 120</code>.</p>
  <p>To do this sort of thing you should have the headers bypass your filtering, which is a pain:</p>
  <pre><code>(samtools view --header-only $file; \
   samtools view --no-header $file | awk …) | samtools view -b -o $outfile
  </code></pre>
  <p>Or add to your awk conditions so that header lines are always kept, which is slightly less painful:</p>
  <pre><code>samtools view --with-header $file | \
  awk '/^@/ || (length($10) &gt;= 120 &amp;&amp; length($10) &lt;= 180)' | \
  samtools view -b -o $outfile
  </code></pre>
  <p>Much easier to use samtools's built-in filtering expressions to do your filtering instead:</p>
  <pre><code>samtools view -e 'qlen &gt;= 120 &amp;&amp; qlen &lt;= 180' -b -o $outfile $file
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-05-24T09:06:33.063Z
Answers.LastActivityDate: 2024-05-26T08:51:12.037Z
Answers.LastEditDate: 2024-05-26T08:51:12.037Z
Answers.UserDisplayName: John Marshall
Body: <p>Hi I'm wondering how I can subset a BAM file based on read length. Precisely I want just read lengths of 120 to 180 bp reads in the new BAM file.</p>
  <p>I'm trying several way, but none of them giving the output</p>
  <pre><code>samtools view -h <span class="math-container">$BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.sort.rmdup.bam | 
  awk '{if (length($10) &gt;= 120 &amp;&amp; length($10) &lt;= 180) print $0}' |
  samtools view -b - &gt; $</span>BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.midSny.subset.bam
  
  [E::sam_parse1] missing SAM header
  [W::sam_read1] parse error at line 1
  [main_samview] truncated file.
  
  samtools view -h <span class="math-container">$BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.sort.rmdup.bam |
  awk 'NR==1 || (length($</span>10) &gt;= 120 &amp;&amp; length(<span class="math-container">$10) &lt;= 180)' |
  samtools view -b - &gt; $</span>BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.midSny.subset.bam
  
  E::sam_parse1] missing SAM header
  [W::sam_read1] parse error at line 2
  [main_samview] truncated file.
  </code></pre>
  <p>However, when I'm providing only single condition then it's working well,</p>
  <pre><code> samtools view -h S10B_ZKDN240002407-1A_H5VMFDSXC_L1_bwaM.sort.rmdup.bam | 
      awk '{if (length(<span class="math-container">$10) &lt;= 180) print $</span>0}' |
      samtools view -b - &gt; temp.bam
  </code></pre>
  <p>But, if I use this temp.bam file for 2nd condition, then it's not working</p>
  <pre><code> samtools view -h temp.bam |
   awk '{if (length(<span class="math-container">$10) &gt;= 120) print $</span>0}'  |
   samtools view -b - &gt;  S10B_ZKDN240002407-1A_H5VMFDSXC_L1_test2.midSny.subset.bam
  
  [E::sam_parse1] missing SAM header
  [W::sam_read1] parse error at line 1
  [main_samview] truncated file. 
  </code></pre>
  
Tags: <sam><genomics><samtools><subset>
Title: How to subset BAM file based on read length range (120-180) bp?
ThreadText: # How to subset BAM file based on read length range (120-180) bp?
  
  (Asked by: username_0 on 2024-05-23)
  
  Hi I'm wondering how I can subset a BAM file based on read length. Precisely I want just read lengths of 120 to 180 bp reads in the new BAM file.
  
  I'm trying several way, but none of them giving the output
  
  ```
  samtools view -h $BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.sort.rmdup.bam | 
  awk '{if (length($10) >= 120 && length($10) <= 180) print $0}' |
  samtools view -b - > $BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.midSny.subset.bam
  
  [E::sam_parse1] missing SAM header
  [W::sam_read1] parse error at line 1
  [main_samview] truncated file.
  
  samtools view -h $BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.sort.rmdup.bam |
  awk 'NR==1 || (length($10) >= 120 && length($10) <= 180)' |
  samtools view -b - > $BAM/S14B_ZKDN240002415-1A_H5VMFDSXC_L1_bwaM.midSny.subset.bam
  
  E::sam_parse1] missing SAM header
  [W::sam_read1] parse error at line 2
  [main_samview] truncated file.
  
  ```
  
  However, when I'm providing only single condition then it's working well,
  
  ```
   samtools view -h S10B_ZKDN240002407-1A_H5VMFDSXC_L1_bwaM.sort.rmdup.bam | 
      awk '{if (length($10) <= 180) print $0}' |
      samtools view -b - > temp.bam
  
  ```
  
  But, if I use this temp.bam file for 2nd condition, then it's not working
  
  ```
   samtools view -h temp.bam |
   awk '{if (length($10) >= 120) print $0}'  |
   samtools view -b - >  S10B_ZKDN240002407-1A_H5VMFDSXC_L1_test2.midSny.subset.bam
  
  [E::sam_parse1] missing SAM header
  [W::sam_read1] parse error at line 1
  [main_samview] truncated file. 
  
  ```
  
  **Tags:** sam, genomics, samtools, subset
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2024-05-24. Score: 3
  
  Your file's headers are also being passed through `awk`. Header lines' fields are all fairly short (and these lines may have fewer than ten fields, so `length($10)` would be zero), so they are all being removed by the awk conditions of the form `length($10) >= 120`.
  
  To do this sort of thing you should have the headers bypass your filtering, which is a pain:
  
  ```
  (samtools view --header-only $file; \
   samtools view --no-header $file | awk …) | samtools view -b -o $outfile
  
  ```
  
  Or add to your awk conditions so that header lines are always kept, which is slightly less painful:
  
  ```
  samtools view --with-header $file | \
  awk '/^@/ || (length($10) >= 120 && length($10) <= 180)' | \
  samtools view -b -o $outfile
  
  ```
  
  Much easier to use samtools's built-in filtering expressions to do your filtering instead:
  
  ```
  samtools view -e 'qlen >= 120 && qlen <= 180' -b -o $outfile $file
  
  ```
  
  
  ---
  
  