<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>

CC BY-SA 4.0
2021-07-16T21:38:26.423Z
2021-07-20T19:07:38.047Z
2021-07-20T19:07:38.047Z
opplatek
<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>

CC BY-SA 4.0
2019-08-19T11:27:11.207Z
The problem with macs2 if you look for cutting sites is that they throw away all the mates (which is 50% of you reads) + only keeps properly aligned. So that is at least 50%. I'll take a look at changing it into bed and see how that works. Thanks
Maarten-vd-Sande
CC BY-SA 4.0
2019-08-19T10:29:23.910Z
2019-08-19T10:29:23.910Z
StupidWolf
<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>

CC BY-SA 4.0
2019-08-19T08:00:17.577Z
That error is probably from Picard, making it irrelevant. Your script should work fine for macs2.
Devon Ryan
CC BY-SA 4.0
2019-08-20T11:46:34.877Z
@DevonRyan Yep, my bad. Everything works just fine as is with macs2. Still weird that picard does not accept the file, but that is another problem.
Maarten-vd-Sande
CC BY-SA 4.0
2019-08-20T13:15:26.060Z
Picard is incredibly picky :(
Devon Ryan
<pysam><macs2>
Convert paired-end BAM into a single-end BAM and keep all the reads
# Convert paired-end BAM into a single-end BAM and keep all the reads

(Asked by: username_0 on 2019-08-19)

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)

**Tags:** pysam, macs2
> **Comments:**
> **username_1**: That error is probably from Picard, making it irrelevant. Your script should work fine for macs2.  
> **username_0**: @DevonRyan Yep, my bad. Everything works just fine as is with macs2. Still weird that picard does not accept the file, but that is another problem.  
> **username_1**: Picard is incredibly picky :(  

---

## 2 answers

### Answer by username_2 on 2021-07-16. Score: 2

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.


---

### Answer by username_3 on 2019-08-19. Score: 2

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  

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.

> **Comments:**
> **username_0**: The problem with macs2 if you look for cutting sites is that they throw away all the mates (which is 50% of you reads) + only keeps properly aligned. So that is at least 50%. I'll take a look at changing it into bed and see how that works. Thanks  

---

