Homec4science
Diffusion scATAC-seq (master) Browse

scATAC-seq/scripts/bam_tools/split_in_two.py
faf6f1d50a68master

Last Change
Raw File

split_in_two.py
View Options

import optparse
import sys
import os
import typing as tp
import pysam
def split_bam(file_in, file_out):
bam_in = pysam.AlignmentFile(file_in)
bam_out = pysam.AlignmentFile(file_out, template=bam_in, mode="wb")
for read in bam_in:
# check whether there is a pair
read_start = read.reference_start
read_flags = list(str(bin(read.flag))[2:][::-1])
read_is_paired = read_flags[0] == '1'
read_is_rev = read_flags[4] == '1'
read_is_first_in_pair = read_flags[6] == '1'
read_name = read.query_name
mate_start = read.next_reference_start
mate_is_rev = read_flags[5] == '1'
# check that read and fragment are OK
# qc
if read.is_qcfail:
continue
# check pair
elif not read_is_paired:
continue
# --> -->
elif(not read_is_rev and not mate_is_rev):
continue
# <-- <--
elif (read_is_rev and read_is_rev):
continue
# <-- -->
elif (read_is_rev and not mate_is_rev) and (read_start < mate_read_start):
continue
# <-- -->
elif (not read_is_rev and mate_is_rev) and (read_start > mate_start):
continue
# Split the fragment in two equally long fragments.
# Each read has a length of 1.
# Reads that do not create a proper fragment
# with their pair read are filtered (can only
# create 2 fragments from a fragment!).
if read_is_first_in_pair:
# strand related parameters
# r1 is fw
# r1 r3
# ---> --->
# |---------|----------|
# <--- <---
# r4 r2
if not read_is_rev:
frag_ref_id = read.reference_id
frag_start = read.reference_start
frag_len = read.template_length
frag_end = frag_start + frag_len - 1
frag_mid = frag_start + (frag_len // 2)
frag_cell = read.get_tag("CB")
# read 1
r1_ref_id = frag_ref_id
r1_next_ref_id = r1_ref_id
r1_start = frag_start # read start inclusive
r1_flags = ['0' for _ in range(12)]
r1_flags[0] = '1' # paired
r1_flags[1] = '1' # proper pair
r1_flags[5] = '1' # mate is rev
r1_flags[6] = '1' # 1st in pair
r1_name = "%s_r1" % read_name
r1_tags = (("CB", frag_cell),)
r1_query_seq = 'N'
# read 2
r2_ref_id = frag_ref_id
r2_next_ref_id = r2_ref_id
r2_start = frag_end -1 # read start inclusive
r2_flags = ['0' for _ in range(12)]
r2_flags[0] = '1' # paired
r2_flags[1] = '1' # proper pair
r2_flags[4] = '1' # read is rev
r2_flags[7] = '1' # mate 1st in pair
r2_name = "%s_r2" % read_name
r2_tags = (("CB", frag_cell),)
r2_query_seq = 'N'
# read 3
r3_ref_id = frag_ref_id
r3_next_ref_id = r3_ref_id
r3_start = frag_mid # read start inclusive
r3_flags = ['0' for _ in range(12)]
r3_flags[0] = '1' # paired
r3_flags[1] = '1' # proper pair
r3_flags[5] = '1' # mate is rev
r3_flags[6] = '1' # 1st in pair
r3_name = "%s_r3" % read_name
r3_tags = (("CB", frag_cell),)
r3_query_seq = 'N'
# read 4
r4_ref_id = frag_ref_id
r4_next_ref_id = r4_ref_id
r4_start = frag_mid - 1 # read start inclusive
r4_flags = ['0' for _ in range(12)]
r4_flags[0] = '1' # paired
r4_flags[1] = '1' # proper pair
r4_flags[4] = '1' # read is rev
r4_flags[7] = '1' # mate 1st in pair
r4_name = "%s_r4" % read_name
r4_tags = (("CB", frag_cell),)
r4_query_seq = 'N'
# fragment lengths
frag14_len = r4_start - r1_start + 1
frag23_len = r2_start - r3_start + 1
r1_tlen = frag14_len
r2_tlen = -frag23_len
r3_tlen = frag23_len
r4_tlen = -frag14_len
# r1 is rv
# r2 r4
# ---> --->
# |---------|----------|
# <--- <---
# r3 r1
else:
frag_ref_id = read.reference_id
frag_start = read.next_reference_start
frag_len = abs(read.template_length)
frag_end = frag_start + frag_len - 1
frag_mid = frag_start + (frag_len // 2)
frag_cell = read.get_tag("CB")
# read 1
r1_ref_id = frag_ref_id
r1_next_ref_id = r1_ref_id
r1_start = frag_end - 1 # read start inclusive
r1_flags = ['0' for _ in range(12)]
r1_flags[0] = '1' # paired
r1_flags[1] = '1' # proper pair
r1_flags[4] = '1' # read is rev
r1_flags[6] = '1' # 1st in pair
r1_name = "%s_r1" % read_name
r1_tags = (("CB", frag_cell),)
r1_query_seq = 'N'
# read 2
r2_ref_id = frag_ref_id
r2_next_ref_id = r2_ref_id
r2_start = frag_start # read start inclusive
r2_flags = ['0' for _ in range(12)]
r2_flags[0] = '1' # paired
r2_flags[1] = '1' # proper pair
r2_flags[5] = '1' # mate is rev
r2_flags[7] = '1' # mate 1st in pair
r2_name = "%s_r2" % read_name
r2_tags = (("CB", frag_cell),)
r2_query_seq = 'N'
# read 3
r3_ref_id = frag_ref_id
r3_next_ref_id = r3_ref_id
r3_start = frag_mid - 1 # read start inclusive
r3_flags = ['0' for _ in range(12)]
r3_flags[0] = '1' # paired
r3_flags[1] = '1' # proper pair
r3_flags[4] = '1' # read is rev
r3_flags[6] = '1' # 1st in pair
r3_name = "%s_r3" % read_name
r3_tags = (("CB", frag_cell),)
r3_query_seq = 'N'
# read 4
r4_ref_id = frag_ref_id
r4_next_ref_id = r4_ref_id
r4_start = frag_mid # read start inclusive
r4_flags = ['0' for _ in range(12)]
r4_flags[0] = '1' # paired
r4_flags[1] = '1' # proper pair
r4_flags[5] = '1' # mate is rev
r4_flags[7] = '1' # mate 1st in pair
r4_name = "%s_r4" % read_name
r4_tags = (("CB", frag_cell),)
r4_query_seq = 'N'
# fragment lengths
frag14_len = r4_start - r1_start + 1
frag23_len = r2_start - r3_start + 1
r1_tlen = -frag14_len
r2_tlen = frag23_len
r3_tlen = -frag23_len
r4_tlen = frag14_len
# create the reads
read1 = pysam.AlignedSegment()
read1.query_name = r1_name
read1.flag = int(''.join(r1_flags)[::-1], base=2)
read1.reference_id = frag_ref_id
read1.reference_start = r1_start
read1.next_reference_id = frag_ref_id
read1.next_reference_start = r4_start
read1.tags = r1_tags
read1.template_length = r1_tlen
read1.query_sequence = r1_query_seq
read2 = pysam.AlignedSegment()
read2.query_name = r2_name
read2.flag = int(''.join(r2_flags)[::-1], base=2)
read2.reference_id = frag_ref_id
read2.reference_start = r2_start
read2.next_reference_id = frag_ref_id
read2.next_reference_start = r3_start
read2.tags = r2_tags
read2.template_length = r2_tlen
read2.query_sequence = r2_query_seq
read3 = pysam.AlignedSegment()
read3.query_name = r3_name
read3.flag = int(''.join(r3_flags)[::-1], base=2)
read3.reference_id = frag_ref_id
read3.reference_start = r3_start
read3.next_reference_id = frag_ref_id
read3.next_reference_start = r2_start
read3.tags = r3_tags
read3.template_length = r3_tlen
read3.query_sequence = r3_query_seq
read4 = pysam.AlignedSegment()
read4.query_name = r4_name
read4.flag = int(''.join(r4_flags)[::-1], base=2)
read4.reference_id = frag_ref_id
read4.reference_start = r4_start
read4.next_reference_id = frag_ref_id
read4.next_reference_start = r1_start
read4.tags = r4_tags
read4.template_length = r4_tlen
read4.query_sequence = r4_query_seq
# write
bam_out.write(read1)
bam_out.write(read2)
bam_out.write(read3)
bam_out.write(read4)
bam_in.close()
bam_out.close()
if __name__ == "__main__":
# parse options
usage = "usage: %s [options]" % os.path.basename(__file__)
epilog = "This program reads a bam file and split the fragments within in two fragments " \
"of equal size." \
"Written by Romain Groux, June 2019"
parser = optparse.OptionParser(usage=usage, epilog=epilog)
parser.add_option("-i", "--input", dest="file_in", default=None, type="string", action="store",
help="the addresse of the bam file to filter.")
parser.add_option("-o", "--output", dest="file_out", default=None, type="string", action="store",
help="The addresse of the output file.")
(options, args) = parser.parse_args()
file_in = options.file_in
file_out = options.file_out
# check options
if file_in is None:
print("Error! No input file given (-i)!", sys.stderr)
exit(1)
elif not os.path.isfile(file_in):
print("Error! %s does not exist!" % file_in)
exit(1)
elif file_out is None:
print("Error! no output file given (-o)!")
split_bam(file_in, file_out)
c4science · Help