#-*-coding:utf-8-*-

import os
import sys
import gzip
import subprocess
import re
import HTSeq
import random
from Bio import SeqIO
from os.path import join

class ReadsAlign:
    def __init__(self, in_fq1 = None, in_fq2 = None, name = None, align_result = None, bam_file = None, out = os.getcwd()):

        self.in_fq1 = in_fq1
        self.in_fq2 = in_fq2
        self.name = name

        self.align_result = align_result

        self.bam = bam_file
        if self.bam is not None:
            so = HTSeq.BAM_Reader(self.bam).get_header_dict().get('HD')['SO']
            if so == 'queryname':
                self.sort_by = 'name'
            elif so == 'coordinate':
                self.sort_by = 'pos'
            else:
                sys.exit('BAM file must be sorted.')

        self.out = out

    def align(self, aligner, aligner_path, align_opts, samtools_path, index, thread = 5, sort_by = 'pos'):
        
        if sort_by == 'name':
            sort_opts = '-n'
        elif sort_by == 'pos':
            sort_opts = ''

        aln = join(aligner_path, aligner)
        samtools = join(samtools_path, 'samtools')
        align_cmd = f'{aln} {align_opts} -p {thread} -x {index} -1 {self.in_fq1} -2 {self.in_fq2} | {samtools} sort {sort_opts} -@ {thread} -o {self.out}/{self.name}.sort.bam'

        print(align_cmd)
        align_ret = subprocess.run(align_cmd, shell = True, check = True, stdout = subprocess.PIPE, stderr = subprocess.PIPE, universal_newlines = True)
        print(align_ret)
        self.align_result = align_ret.stderr

        self.bam = '{}/{}.sort.bam'.format(self.out, self.name)
        self.sort_by = sort_by

    def align_stat(self):

        stat_fh = open(join(self.out, self.name + '.map_stat.xls'), 'w')
        stat_fh.write('Sample\tTotal paired reads\tPaired mapped reads\tUnpaired mapped reads\tUnmapped reads\tTotal mapped\n')

        pair_reads = re.search('(\d+)\sreads; of these:', self.align_result, re.M)
        pair_mapped1 = re.search('(\d+)\s.*aligned concordantly exactly 1 time', self.align_result, re.M)
        pair_mapped2 = re.search('(\d+)\s.*aligned concordantly >1 times', self.align_result, re.M)
        pair_mapped3 = re.search('(\d+)\s.*aligned discordantly 1 time', self.align_result, re.M)
        unmapped = re.search('(\d+)\s.*aligned 0 times$', self.align_result, re.M)
        unpair_mapped1 = re.search('(\d+)\s.*aligned exactly 1 time', self.align_result, re.M)
        unpair_mapped2 = re.search('(\d+)\s.*aligned >1 times', self.align_result, re.M)

        if all([pair Limit your response to 100 words.
user