#!/usr/bin/python

"""
   Try to identify G>T and C>T errors in a bamfile.
   This type of error is associated with a library
   preparation problem where the oxidative stress is higher.

   Usage: deepsnvres.dbsnp.tsv path-to-bamfiles
"""

import sys
import pysam

def processBAM(sample, chrx, pos, alt, path):
	bamfile = pysam.AlignmentFile(path + "/" + sample + ".bam", "rb")
	read1 = {"A": 0, "C": 0, "T": 0, "G": 0, "N": 0}
	read2 = {"A": 0, "C": 0, "T": 0, "G": 0, "N": 0}
	for read in bamfile.fetch(chrx, pos, pos + 1):
		seq    = read.query_sequence
		alipos = read.get_aligned_pairs()
		for i in alipos:
			if i[1] == pos and i[0] != None:
				if read.is_read1:
					read1[seq[i[0]]] += 1
				if read.is_read2:
					read2[seq[i[0]]] += 1
	bamfile.close()
	print "\t".join([sample, chrx, str(pos), alt, str(read1["A"]), str(read1["C"]), str(read1["G"]), str(read1["T"]), str(read2["A"]), str(read2["C"]), str(read2["G"]), str(read2["T"])])

pathtobam = sys.argv[2]
print "sample\tchr\tpos\talt\tread1A\tread1C\tread1G\tread1T\tread2A\tread2C\tread2G\tread2T"

for i in open(sys.argv[1]): # read deepsnvres.tsv
	fields = i.rstrip().split("\t")
	sample = fields[0]
	chrx   = fields[1]
	pos    = fields[2]
	ref    = fields[3]
	alt    = fields[4]

	if (ref == "C" and alt == "A") or (ref == "G" and alt == "T"):
		# this is what are we looking for
		pos = int(pos) - 1 # deepsnvres 1-based, pysam 0-based
		processBAM(sample, chrx, pos, alt, pathtobam)
