text: user5084: I want to call peaks using MACS2 on my paired-end ATAC-seq BAM and I am interested in cutting sites. So I make use of the -BAM option (opposed to the -BAMPE) option. However this throws away all the mates, and all the reads of which the mate did not align (that's more than 50%!!!). My idea was to align paired-end, and then remove all the information from the BAM that it is paired-end. So my question is how do I:
  Convert paired-end BAM into a single-end BAM and keep all the reads
  I made a start based on the SAM-format specification , but when I run the result of this script through markduplicates I get the error:
  ERROR: Record 1, Read name SRR6245286.60336908., MRNM should not be set for unpaired read.
  
  I can not find what this MRNM is...
  This is my code I used to generate the faulty 'single-end' bam, but I am open to suggestions (e.g. bash, preferably no R). 
  import pysam
  
  paired =         1
  proper_pair =    2
  mate_unmapped =  8
  mate_reverse =   32
  first_in_pair =  64
  second_in_pair = 128
  
  bam_in = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7.samtools-coordinate.bam", "rb")
  bam_out = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7_SE.samtools-coordinate.bam", "wb", template=bam_in)
  for line in bam_in:
      if line.flag & second_in_pair:
          line.pos += line.template_length
      line.next_reference_id = 0
      line.next_reference_start = 0
      line.flag &= ~(paired + proper_pair + mate_unmapped + mate_reverse + first_in_pair + second_in_pair)
  
      bam_out.write(line)
  
  
  user5077: Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly. 
  You can also use samtools rmdup too
  samtools rmdup <bamfile> <output>
  
  Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.
  I would suggest removing the duplicates and then converting the bam file to a bed file. Macs2 can take in the bed format. Any extremely crude python script that converts a bam into a bed:
  import pysam
  samfile = pysam.AlignmentFile("test.bam", "rb")
  open("test.bed", 'w').close()
  f = open("test.bed", "a")
  STRAND = ["+","-"]
  for read in samfile.fetch():
      STR = STRAND[int(read.is_reverse)]
      BED =[read.reference_name,read.pos,read.reference_end,".",read.mapq,STR]
      f.write('\t'.join([str(i) for i in BED]))
      f.write('\n')
  f.close()
  
  You can also try the conversion with bedtools or definitely refine that script above with pandas etc. 
  
  user486: MRNM stands for "Mate reference index". So Picard found something in the RNEXT field which should be set only for paired-end reads but the rest of the file looks like single-end.
  The problematic line in your code is:
  line.next_reference_id = 0
  
  This sets the RNEXT SAM field to whatever Pysam stores as a reference with index 0 (next_reference_id). This is not the same as reference name (next_reference_name). In single-end SAM file this field should be set to *. You cannot set * with next_reference_id but you can set it with next_reference_name:
  line.next_reference_name = "*" 
  
  I would also recommend to change next_reference_start line to:
  line.next_reference_start = -1 
  
  just to be sure Pysam handles this field as not set (I tried to put there 0 Pysam changed it to 1 for some reason). I would also add line:
  line.template_length = -1
  
  because this field should also be empty in single-end SAM.
  This worked in my case. If Picard keeps complaining you can add VALIDATION_STRINGENCY=LENIENT or SILENT to your Picard command to suppress exiting of the command at errors.
  
original_text: 5084: <p>I want to call peaks using MACS2 on my paired-end ATAC-seq BAM and I am interested in cutting sites. So I make use of the -BAM option (opposed to the -BAMPE) option. However this throws away all the mates, and all the reads of which the mate did not align (that's more than 50%!!!). My idea was to align paired-end, and then remove all the information from the BAM that it is paired-end. So my question is how do I:</p>
  
  <p><strong>Convert paired-end BAM into a single-end BAM and keep all the reads</strong></p>
  
  <p>I made a start based on the <a href="https://samtools.github.io/hts-specs/SAMv1.pdf" rel="nofollow noreferrer">SAM-format specification </a>, but when I run the result of this script through markduplicates I get the error:</p>
  
  <pre><code>ERROR: Record 1, Read name SRR6245286.60336908., MRNM should not be set for unpaired read.
  </code></pre>
  
  <p>I can not find what this MRNM is...</p>
  
  <p>This is my code I used to generate the faulty <em>'single-end'</em> bam, but I am open to suggestions (e.g. bash, preferably no R). </p>
  
  <pre><code>import pysam
  
  paired =         1
  proper_pair =    2
  mate_unmapped =  8
  mate_reverse =   32
  first_in_pair =  64
  second_in_pair = 128
  
  bam_in = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7.samtools-coordinate.bam", "rb")
  bam_out = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7_SE.samtools-coordinate.bam", "wb", template=bam_in)
  for line in bam_in:
      if line.flag &amp; second_in_pair:
          line.pos += line.template_length
      line.next_reference_id = 0
      line.next_reference_start = 0
      line.flag &amp;= ~(paired + proper_pair + mate_unmapped + mate_reverse + first_in_pair + second_in_pair)
  
      bam_out.write(line)
  </code></pre>
   5077: <p>Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly. </p>
  
  <p>You can also use samtools rmdup too</p>
  
  <pre><code>samtools rmdup &lt;bamfile&gt; &lt;output&gt;
  </code></pre>
  
  <p>Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.</p>
  
  <p>I would suggest removing the duplicates and then converting the bam file to a bed file. Macs2 can take in the bed format. Any extremely crude python script that converts a bam into a bed:</p>
  
  <pre><code>import pysam
  samfile = pysam.AlignmentFile("test.bam", "rb")
  open("test.bed", 'w').close()
  f = open("test.bed", "a")
  STRAND = ["+","-"]
  for read in samfile.fetch():
      STR = STRAND[int(read.is_reverse)]
      BED =[read.reference_name,read.pos,read.reference_end,".",read.mapq,STR]
      f.write('\t'.join([str(i) for i in BED]))
      f.write('\n')
  f.close()
  </code></pre>
  
  <p>You can also try the conversion with bedtools or definitely refine that script above with pandas etc. </p>
   486: <p>MRNM stands for &quot;Mate reference index&quot;. So Picard found something in the RNEXT field which should be set only for paired-end reads but the rest of the file looks like single-end.</p>
  <p>The problematic line in your code is:</p>
  <pre><code>line.next_reference_id = 0
  </code></pre>
  <p>This sets the RNEXT SAM field to whatever Pysam stores as a reference with index 0 (<code>next_reference_id</code>). This is <strong>not</strong> the same as reference name (<code>next_reference_name</code>). In single-end SAM file this field should be set to <code>*</code>. You cannot set <code>*</code> with <code>next_reference_id</code> but you can set it with <code>next_reference_name</code>:</p>
  <pre><code>line.next_reference_name = &quot;*&quot; 
  </code></pre>
  <p>I would also recommend to change <code>next_reference_start</code> line to:</p>
  <pre><code>line.next_reference_start = -1 
  </code></pre>
  <p>just to be sure Pysam handles this field as not set (I tried to put there <code>0</code> Pysam changed it to <code>1</code> for some reason). I would also add line:</p>
  <pre><code>line.template_length = -1
  </code></pre>
  <p>because this field should also be empty in single-end SAM.</p>
  <p>This worked in my case. If Picard keeps complaining you can add <code>VALIDATION_STRINGENCY=LENIENT</code> or <code>SILENT</code> to your Picard command to suppress exiting of the command at errors.</p>
  