Answers.Body: <p>For at least <strong>one in a pair</strong> matching the sequence, if the read names of pairs are identical, some grep magic could do the job:</p>
  
  <ol>
  <li>Print the sequence at the beginning of the line and read the name in the second column</li>
  <li>Find all lines starting with the sequence</li>
  <li>Extract only read names</li>
  </ol>
  
  <p>And then just grep all found reads and save into the BAM file:</p>
  
  <pre><code>SEQ=&lt;my_sequence&gt;
  samtools view unfiltered.bam | \
     awk '{print $10 "\t" $1}' | \
     grep "^$SEQ" | \
     cut -f 2 &gt; list_of_matching_names
  samtools view unfiltered.bam | \
     grep -f list_of_matching_names | \
     samtools view -b &gt; filtered.bam
  </code></pre>
  
  <p>Two small notes:</p>
  
  <ol>
  <li>You might want to consider also reverse complementary sequences.</li>
  <li>You also might want to keep header of BAM file (you can extract it at the beginning and add it in the end).</li>
  </ol>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-10T14:03:42.830Z
Answers.LastActivityDate: 2017-11-10T16:53:20.583Z
Answers.LastEditDate: 2017-11-10T16:53:20.583Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <p><strong>First option</strong></p>
  
  <p>My solution uses the following steps:</p>
  
  <ul>
  <li>use picard sortsam to sort the records on query-name (not samtools sort because the order <a href="https://github.com/samtools/hts-specs/issues/5" rel="nofollow noreferrer">is not the same</a> between java and C ) </li>
  <li>use <a href="https://docs.oracle.com/javase/8/docs/technotes/tools/unix/jjs.html" rel="nofollow noreferrer">jjs</a> (java scripting engine) and the <a href="https://github.com/samtools/htsjdk" rel="nofollow noreferrer">htsjdk library</a> to build a <code>buffer</code>of reads having the same name. If any read starts with a pattern, print the whole <code>buffer</code>. The htsjdk should be already contained in the picard.jar.</li>
  <li>sort back to coordinate order using samtools.</li>
  </ul>
  
  <p>All in one:</p>
  
  <pre><code>java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
  jjs -cp /path/to/picard.jar  script.js |\
  samtools sort -T TMP -o out.bam - 
  </code></pre>
  
  <p>with script.js</p>
  
  <pre><code>var out=java.lang.System.out;
  var File = Java.type("java.io.File");
  var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
  var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
  var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
  var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
  var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
  var iter=samReader.iterator();
  var buffer=[];
  for(;;)
          {
          var rec=null;
          if(iter.hasNext()) rec=iter.next();
          if(rec==null || (buffer.length&gt;0 &amp;&amp; !rec.getReadName().equals(buffer[0].getReadName())))
              {
              var i=0;
              for(i=0;i&lt; buffer.length;++i)
                  {
                  if(buffer[i].getReadString().match(/^AAAA/)) break;
                  }
              if(i!=buffer.length)
                  {
                  for(i=0;i&lt; buffer.length ;++i)
                      {
                      samWriter.addAlignment(buffer[i]);
                      }
                  }
              if(rec==null) break;
              buffer=[];
              }
          buffer.push(rec);
          }
  
  samReader.close();
  samWriter.close();
  </code></pre>
  
  <p><strong>A second option</strong> is to use samtools view and awk to export the read names:</p>
  
  <pre><code> samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq &gt; names.txt
  </code></pre>
  
  <p>and to use <a href="http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads" rel="nofollow noreferrer">Picard FilterSamReads</a> and READ_LIST_FILE to filter the reads.</p>
  
  <p><strong>UPDATE: 3rd solution:</strong></p>
  
  <p>I've updated my tool <a href="http://lindenb.github.io/jvarkit/SamJdk.html" rel="nofollow noreferrer">samjdk</a>: there is now a new option <code>--pair</code> to work with query-name-sorted BAM files.</p>
  
  <p>The command would be:</p>
  
  <pre><code>java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
  java -jar dist/samjdk.jar --samoutputformat BAM --pair \
       -e 'return records.stream().anyMatch(R-&gt;R.getReadString().startsWith("AAAAAAG"));' |\
  samtools sort -T TMP -o out.bam - 
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-15T11:08:05.287Z
Answers.Text: Does this account for the fact that BAM stores the sequence of a read in genome orientation, rather than read orientation? So if the sequence of a read mapped to the plus strand starts with AAAA it will be recorded in the BAM as ending with TTTT.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-15T11:17:51.420Z
Answers.Text: @ia no but it would be easy to test this: `(R.getReadNegativeStrandFlag()?R.getReadString().endsWith("TTTT"):R.getReadString().startsWith("AAAA"))`
Answers.UserDisplayName: Pierre
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-15T16:38:50.737Z
Answers.Text: I got it to work now with samjdk as explained in 3rd solution. Thanks a lot!
Answers.UserDisplayName: 719016
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-10T13:54:58.450Z
Answers.LastActivityDate: 2017-11-10T16:32:35.243Z
Answers.LastEditDate: 2017-11-10T16:32:35.243Z
Answers.UserDisplayName: Pierre
Answers.Body: <p>I am one of the authors of the <code>UMI-tools</code> pacakge, and although you probably could hack together a way to do this with our package, its not what its designed for and the solutions above are probably better.</p>
  
  <p>The best way to do this would be to filter the fastq before mapping if thats a possibility, if it has to be the BAM, then here is a python implementation similar to Pierres java solution above. This one fiddles about with strands and cigar operations a bit more. </p>
  
  <p>First sort the reads by name with <code>samtools sort -n</code></p>
  
  <p>then</p>
  
  <pre class="lang-python prettyprint-override"><code>import pysam
  inbam = pysam.AlignmentFile("mybam.bam")
  outbam = pysam.AlignmentFile("mybam.filtered.bam", "wb", template=inbam)
  seq_to_find = "AAGG"
  rc_seq_to_find = "CCTT"
  
  def test_and_out_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find):
      '''This function takes a list of reads and check if any of them
      have the prespecific sequences at their start. If the sequence
      is detected the whole buffer is output. The start of the read 
      is the end of the alignment if the read is on the '-' strand. 
      Also checks that the read doesn't start with hard-clipped sequence'''
  
      do_output=False
  
      for segment in read_buffer:
  
          if segment.is_reverse:
              # the last operation is stored in the first part of the
              #last of the cigar tuples. Operation 5 is Hard-clip
              output_condition = (not segment.cigartuples[-1][0] == 5) and \
                segment.query_sequence.endswith(rc_seq_to_find)
          else:
              output_condition = (not segment.cigartuples[0][0] == 5) and \
                segment.query_sequence.startswith(seq_to_find)
  
          if output_condition:
              do_output=True
              break
  
      if do_output:
         for segment in read_buffer:
             outbam.write(segment)
  
  last_name=""
  read_buffer=[]
  
  # Loop over reads in the file
  for read in inbam.fetch(until_eof=True):
  
      # if the name hasn't changed, the name is part of the 
      # same "pair"
      if read.query_name==last_name:
          buffer.append(read)
          continue
  
      # if we got here then the a new pair is being examined
      # test the old pair and output if neccesary. 
      test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)  
  
      # make a new buffer with just the single new read.
      last_name=read.query_name
      read_buffer=[read]
  
  # the last set of reds will not have been output.  
  test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-15T12:47:09.740Z
Answers.LastActivityDate: 2017-11-15T12:47:09.740Z
Answers.UserDisplayName: Ian Sudbery
Body: <p>What's the most straightforward way to filter a BAM file for read pairs where one or both of the reads starts with a given sequence pattern? Would a tool that deals with UMIs work for this? Which one would be best?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-15T11:13:34.367Z
Comments.Text: Are you interested in READs that start with a particular sequence or ALIGNMENTS that start with a particular sequence, because they are not the same thing for 2 reasons:
  
  1) The sam format records sequences in genome orientation. Thus if the sequence of the read starts AAAA, but the read maps to the '-' strand the BAM seq field will not start AAAA, but end TTTT
  
  2) Mappers may clip or otherwise only map part of a read. When they do this, they might only record part of the sequence of the read in the SEQ field. Thus a line with the CIGAR 5H90M will not contain the first 5nt of the read seq
Comments.UserDisplayName: Ian Sudbery
Tags: <sam><sequence-analysis><filtering>
Title: Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern
ThreadText: # Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern
  
  (Asked by: username_0 on 2017-11-10)
  
  What's the most straightforward way to filter a BAM file for read pairs where one or both of the reads starts with a given sequence pattern? Would a tool that deals with UMIs work for this? Which one would be best?
  
  **Tags:** sam, sequence-analysis, filtering
  
  > **Comments:**
  > **username_1**: Are you interested in READs that start with a particular sequence or ALIGNMENTS that start with a particular sequence, because they are not the same thing for 2 reasons:  1) The sam format records sequences in genome orientation. Thus if the sequence of the read starts AAAA, but the read maps to the '-' strand the BAM seq field will not start AAAA, but end TTTT  2) Mappers may clip or otherwise only map part of a read. When they do this, they might only record part of the sequence of the read in the SEQ field. Thus a line with the CIGAR 5H90M will not contain the first 5nt of the read seq  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_2 on 2017-11-10. Score: 3
  
  For at least **one in a pair** matching the sequence, if the read names of pairs are identical, some grep magic could do the job:
  
  1. Print the sequence at the beginning of the line and read the name in the second column
  2. Find all lines starting with the sequence
  3. Extract only read names
  
  And then just grep all found reads and save into the BAM file:
  
  ```
  SEQ=<my_sequence>
  samtools view unfiltered.bam | \
     awk '{print $10 "\t" $1}' | \
     grep "^$SEQ" | \
     cut -f 2 > list_of_matching_names
  samtools view unfiltered.bam | \
     grep -f list_of_matching_names | \
     samtools view -b > filtered.bam
  
  ```
  
  Two small notes:
  
  1. You might want to consider also reverse complementary sequences.
  2. You also might want to keep header of BAM file (you can extract it at the beginning and add it in the end).
  
  
  ---
  
  ### Answer by username_3 on 2017-11-10. Score: 3
  
  **First option**
  
  My solution uses the following steps:
  
  * use picard sortsam to sort the records on query-name (not samtools sort because the order [is not the same](https://github.com/samtools/hts-specs/issues/5) between java and C )
  * use [jjs](https://docs.oracle.com/javase/8/docs/technotes/tools/unix/jjs.html) (java scripting engine) and the [htsjdk library](https://github.com/samtools/htsjdk) to build a `buffer`of reads having the same name. If any read starts with a pattern, print the whole `buffer`. The htsjdk should be already contained in the picard.jar.
  * sort back to coordinate order using samtools.
  
  All in one:
  
  ```
  java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
  jjs -cp /path/to/picard.jar  script.js |\
  samtools sort -T TMP -o out.bam - 
  
  ```
  
  with script.js
  
  ```
  var out=java.lang.System.out;
  var File = Java.type("java.io.File");
  var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
  var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
  var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
  var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
  var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
  var iter=samReader.iterator();
  var buffer=[];
  for(;;)
          {
          var rec=null;
          if(iter.hasNext()) rec=iter.next();
          if(rec==null || (buffer.length>0 && !rec.getReadName().equals(buffer[0].getReadName())))
              {
              var i=0;
              for(i=0;i< buffer.length;++i)
                  {
                  if(buffer[i].getReadString().match(/^AAAA/)) break;
                  }
              if(i!=buffer.length)
                  {
                  for(i=0;i< buffer.length ;++i)
                      {
                      samWriter.addAlignment(buffer[i]);
                      }
                  }
              if(rec==null) break;
              buffer=[];
              }
          buffer.push(rec);
          }
  
  samReader.close();
  samWriter.close();
  
  ```
  
  **A second option** is to use samtools view and awk to export the read names:
  
  ```
   samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq > names.txt
  
  ```
  
  and to use [Picard FilterSamReads](http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads) and READ\_LIST\_FILE to filter the reads.
  
  **UPDATE: 3rd solution:**
  
  I've updated my tool [samjdk](http://lindenb.github.io/jvarkit/SamJdk.html): there is now a new option `--pair` to work with query-name-sorted BAM files.
  
  The command would be:
  
  ```
  java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
  java -jar dist/samjdk.jar --samoutputformat BAM --pair \
       -e 'return records.stream().anyMatch(R->R.getReadString().startsWith("AAAAAAG"));' |\
  samtools sort -T TMP -o out.bam - 
  
  ```
  
  > **Comments:**
  > **username_1**: Does this account for the fact that BAM stores the sequence of a read in genome orientation, rather than read orientation? So if the sequence of a read mapped to the plus strand starts with AAAA it will be recorded in the BAM as ending with TTTT.  
  > **username_3**: @ia no but it would be easy to test this: `(R.getReadNegativeStrandFlag()?R.getReadString().endsWith("TTTT"):R.getReadString().startsWith("AAAA"))`  
  > **username_0**: I got it to work now with samjdk as explained in 3rd solution. Thanks a lot!  
  
  ---
  
  ### Answer by username_1 on 2017-11-15. Score: 0
  
  I am one of the authors of the `UMI-tools` pacakge, and although you probably could hack together a way to do this with our package, its not what its designed for and the solutions above are probably better.
  
  The best way to do this would be to filter the fastq before mapping if thats a possibility, if it has to be the BAM, then here is a python implementation similar to Pierres java solution above. This one fiddles about with strands and cigar operations a bit more.
  
  First sort the reads by name with `samtools sort -n`
  
  then
  
  ```python
  import pysam
  inbam = pysam.AlignmentFile("mybam.bam")
  outbam = pysam.AlignmentFile("mybam.filtered.bam", "wb", template=inbam)
  seq_to_find = "AAGG"
  rc_seq_to_find = "CCTT"
  
  def test_and_out_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find):
      '''This function takes a list of reads and check if any of them
      have the prespecific sequences at their start. If the sequence
      is detected the whole buffer is output. The start of the read 
      is the end of the alignment if the read is on the '-' strand. 
      Also checks that the read doesn't start with hard-clipped sequence'''
  
      do_output=False
  
      for segment in read_buffer:
  
          if segment.is_reverse:
              # the last operation is stored in the first part of the
              #last of the cigar tuples. Operation 5 is Hard-clip
              output_condition = (not segment.cigartuples[-1][0] == 5) and \
                segment.query_sequence.endswith(rc_seq_to_find)
          else:
              output_condition = (not segment.cigartuples[0][0] == 5) and \
                segment.query_sequence.startswith(seq_to_find)
  
          if output_condition:
              do_output=True
              break
  
      if do_output:
         for segment in read_buffer:
             outbam.write(segment)
  
  last_name=""
  read_buffer=[]
  
  # Loop over reads in the file
  for read in inbam.fetch(until_eof=True):
  
      # if the name hasn't changed, the name is part of the 
      # same "pair"
      if read.query_name==last_name:
          buffer.append(read)
          continue
  
      # if we got here then the a new pair is being examined
      # test the old pair and output if neccesary. 
      test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)  
  
      # make a new buffer with just the single new read.
      last_name=read.query_name
      read_buffer=[read]
  
  # the last set of reds will not have been output.  
  test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)
  
  ```
  
  
  ---
  
  