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)