1
$\begingroup$

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)
$\endgroup$
3
  • 1
    $\begingroup$ That error is probably from Picard, making it irrelevant. Your script should work fine for macs2. $\endgroup$
    – Devon Ryan
    Aug 19 '19 at 8:00
  • $\begingroup$ @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. $\endgroup$ Aug 20 '19 at 11:46
  • $\begingroup$ Picard is incredibly picky :( $\endgroup$
    – Devon Ryan
    Aug 20 '19 at 13:15
2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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 $\endgroup$ Aug 19 '19 at 11:27
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.