#!/usr/bin/python

import os
import re
import platform

from Bio import AlignIO
from Bio.Seq import Seq
from Bio.Align import AlignInfo
from Bio.Align.Applications import MuscleCommandline

# define constants
HEADER = 1
SEQUENCE = 2

FORWARD = 1
REVERSE = 2

FA_LINES_PER_READ = 2
FQ_LINES_PER_READ = 4
class SequenceData:
    def __init__(self, read_type, telomeric_repeat, interleaved_read_filepath=None, r1_read_filepath=None, r2_read_filepath=None, out_directory=None):
        self._read_type = read_type

        self._interleaved_reads_filepath = interleaved_read_filepath 
        self._r1_reads_filepath = r1_read_filepath 
        self._r2_reads_filepath = r2_read_filepath

        self._out_directory = out_directory
        
        if self._out_directory == None:
            self._out_directory = os.path.abspath(os.path.join(os.path.expanduser('~'), 'terminus_output'))
            os.makedirs(os.path.dirname(self._out_directory), exist_ok=True)

        self._3prime_telomeric_repeat_regex = self._generate_regular_expression(telomeric_repeat, "3'")
        self._5prime_telomeric_repeat_regex = self._generate_regular_expression(str(Seq(telomeric_repeat).reverse_complement()), "5'")

        self._r1_telomeric_reads_filepath = None
        self._r2_telomeric_reads_filepath = None

        self._r1_telomeric_reads_clusters_info_filepath = None
        self._r2_telomeric_reads_clusters_info_filepath = None

        self._r1_telomeric_reads_clusters_directory = None
        self._r2_telomeric_reads_clusters_directory = None
        
        self._r1_telomeric_reads_clusters_alignments_directory = None
        self._r2_telomeric_reads_clusters_alignments_directory = None

        self._r1_telomeric_reads_clusters_consensus_filepath = None
        self._r2_telomeric_reads_clusters_consensus_filepath = None

        self._r1_subtelomeric_reads_filepath = None
        self._r2_subtelomeric_reads_filepath = None

        self._r1_subtelomeric_reads_clusters_info_filepath = None
        self._r2_subtelomeric_reads_clusters_info_filepath = None

        self._r1_subtelomeric_reads_clusters_directory = None
        self._r2_subtelomeric_reads_clusters_directory = None

        self._r1_subtelomeric_reads_clusters_alignments_directory = None
        self._r2_subtelomeric_reads_clusters_alignments_directory = None

        self._r1_subtelomeric_reads_clusters_consensus_filepath = None
        self._r2_subtelomeric_reads_clusters_consensus_filepath = None
    
    def get_read_type(self):
        return self._read_type

    def get_interleaved_reads_filepath(self):
        return self._interleaved_reads_filepath

    def get_r1_reads_filepath(self):
        return self._r1_reads_filepath

    def get_r2_reads_filepath(self):
        return self._r2_reads_filepath
    
    def get_out_directory(self):
        return self._out_directory
    
    def get_3prime_telomeric_repeat_regex(self):
        return self._3prime_telomeric_repeat_regex

    def get_5prime_telomeric_repeat_regex(self):
        return self._5prime_telomeric_repeat_regex
    
    def get_r1_telomeric_reads_filepath(self):
        return self._r1_telomeric_reads_filepath
    
    def get_r2_telomeric_reads_filepath(self):
        return self._r2_telomeric_reads_filepath
    
    def get_r1_telomeric_reads_clusters_info_filepath(self):
        return self._r1_telomeric_reads_clusters_info_filepath
    
    def get_r2_telomeric_reads_clusters_info_filepath(self):
        return self._r2_telomeric_reads_clusters_info_filepath

    def get_r1_telomeric_reads_clusters_directory(self):
        return self._r1_telomeric_reads_clusters_directory
    
    def get_r2_telomeric_reads_clusters_directory(self):
        return self._r2_telomeric_reads_clusters_directory

    def get_r1_telomeric_reads_clusters_alignments_directory(self):    
        return self._r1_telomeric_reads_clusters_alignments_directory
        
    def get_r2_telomeric_reads_clusters_alignments_directory(self):    
        return self._r2_telomeric_reads_clusters_alignments_directory

    def get_r1_telomeric_reads_clusters_consensus_filepath(self): 
        return self._r1_telomeric_reads_clusters_consensus_filepath
        
    def get_r2_telomeric_reads_clusters_consensus_filepath(self):
        return self._r2_telomeric_reads_clusters_consensus_filepath
    
    def get_r1_subtelomeric_reads_filepath(self):
        return self._r1_subtelomeric_reads_filepath
    
    def get_r2_subtelomeric_reads_filepath(self):
        return self._r2_subtelomeric_reads_filepath
    
    def get_r1_subtelomeric_reads_clusters_info_filepath(self):
        return self._r1_subtelomeric_reads_clusters_info_filepath
    
    def get_r2_subtelomeric_reads_clusters_info_filepath(self):
        return self._r2_subtelomeric_reads_clusters_info_filepath

    def get_r1_subtelomeric_reads_clusters_directory(self):
        return self._r1_subtelomeric_reads_clusters_directory
    
    def get_r2_subtelomeric_reads_clusters_directory(self):
        return self._r2_subtelomeric_reads_clusters_directory

    def get_r1_subtelomeric_reads_clusters_alignments_directory(self):    
        return self._r1_subtelomeric_reads_clusters_alignments_directory
        
    def get_r2_subtelomeric_reads_clusters_alignments_directory(self):    
        return self._r2_subtelomeric_reads_clusters_alignments_directory

    def get_r1_subtelomeric_reads_clusters_consensus_filepath(self): 
        return self._r1_subtelomeric_reads_clusters_consensus_filepath
        
    def get_r2_subtelomeric_reads_clusters_consensus_filepath(self):
        return self._r2_subtelomeric_reads_clusters_consensus_filepath
    
    def extract_telomeric_reads(self):
        telomeric_reads_directory = os.path.abspath(os.path.join(self._out_directory, 'telomeric_reads'))
        os.makedirs(telomeric_reads_directory, exist_ok=True)

        self._r1_telomeric_reads_filepath = os.path.abspath(os.path.join(telomeric_reads_directory, 'r1_telomeric_reads.fasta'))
        r1_telomeric_reads_output = open(self._r1_telomeric_reads_filepath, "w+")

        self._r2_telomeric_reads_filepath = os.path.abspath(os.path.join(telomeric_reads_directory, 'r2_telomeric_reads.fasta'))
        r2_telomeric_reads_output = open(self._r2_telomeric_reads_filepath, "w+")

        if self._read_type == "interleaved":
            print("Extracting telomeric reads from an interleaved FASTQ dataset...")
            self._process_reads(self._interleaved_reads_filepath, r1_telomeric_reads_output, r2_telomeric_reads_output)        
            print("Done.\n")

        if self._read_type == "seperated":
            print("Extracting telomeric reads from R1 FASTQ dataset...")
            self._process_reads(input_filepath=self._r1_reads_filepath, r1_output=r1_telomeric_reads_output, r2_output=None)
            print("Done.\n")

            print("Extracting telomeric reads from R2 FASTQ dataset...")
            self._process_reads(input_filepath=self._r2_reads_filepath, r1_output=None, r2_output=r2_telomeric_reads_output)
            print("Done.\n")       

        print("R1 telomeric reads saved to " + self._r1_telomeric_reads_filepath)
        print("R2 telomeric reads saved to " + self._r2_telomeric_reads_filepath)
        print()
 
    def cluster_telomeric_reads(self):   
        # install and make wcd for clustering
        print("Installing and preparing wcdest for clustering...")
        os.chdir("..")
        os.system("git submodule update --init --recursive >/dev/null 2>&1")
        os.chdir("./lib/wcdest/code/")
        os.system("./configure >/dev/null 2>&1")
        # NOTE: texinfo/makeinfo must be installed for wcdest to compile correctly!! 
        os.system("make clean -C ./ >/dev/null 2>&1")
        os.system("make -C ./ >/dev/null 2>&1")
        os.system("make install -C ./ >/dev/null 2>&1")
        os.chdir("../../../src")
        print("Done.\n")

        telomeric_reads_cluster_info_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'cluster_info', 'telomeric_reads'))
        os.makedirs(telomeric_reads_cluster_info_directory, exist_ok=True)
        
        self._r1_telomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(telomeric_reads_cluster_info_directory, 'r1_telomeric_reads_clusters_info.txt'))
        self._r2_telomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(telomeric_reads_cluster_info_directory, 'r2_telomeric_reads_clusters_info.txt'))

        # cluster forward telomeric reads and write results to file 
        print("Clustering telomeric reads in '" + self._r1_telomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r1_telomeric_reads_clusters_info_filepath + " -c " + self._r1_telomeric_reads_filepath + " >/dev/null 2>&1")
        self._r1_telomeric_reads_clusters_directory = self._generate_clustered_reads(self._r1_telomeric_reads_filepath, self._r1_telomeric_reads_clusters_info_filepath, "telomeric", "forward")
        print("Done.\n")

        # cluster reverse telomeric reads and write results to file 
        print("Clustering telomeric reads in '" + self._r2_telomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r2_telomeric_reads_clusters_info_filepath + " -c " + self._r2_telomeric_reads_filepath + " >/dev/null 2>&1")
        self._r2_telomeric_reads_clusters_directory = self._generate_clustered_reads(self._r2_telomeric_reads_filepath, self._r2_telomeric_reads_clusters_info_filepath, "telomeric", "reverse")
        print("Done.\n")
        
        print("R1 telomeric reads clustering results saved to '" + self._r1_telomeric_reads_clusters_directory)
        print("R2 telomeric reads clustering results saved to '" + self._r2_telomeric_reads_clusters_directory)
        print()

    def extract_subtelomeric_reads(self):
        # given forward telomeric clusters, retrieve the subtelomeric sequence reads
        #  from the reverse read file  
        print("Extracting subtelomeric reads from R2 sequence reads file for clustered R1 telomeric reads...")
        for reads_file in os.listdir(self._r1_telomeric_reads_clusters_directory):
            if reads_file.endswith(".fasta"):
                if self._read_type == "interleaved":
                    self._r2_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._interleaved_reads_filepath, os.path.abspath(os.path.join(self._r1_telomeric_reads_clusters_directory, reads_file)), "reverse")
                if self._read_type == "seperated":
                    self._r2_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._r2_reads_filepath, os.path.abspath(os.path.join(self._r1_telomeric_reads_clusters_directory, reads_file)), "reverse")
        print("Done.\n")

        # given reverse telomeric clusters, retrieve the subtelomeric sequence reads
        #  from the forward read file  
        print("Extracting subtelomeric reads from R1 sequence reads file for clustered R2 telomeric reads...")
        for reads_file in os.listdir(self._r2_telomeric_reads_clusters_directory):
            if reads_file.endswith(".fasta"): 
                if self._read_type == "interleaved":
                    self._r1_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._interleaved_reads_filepath, os.path.abspath(os.path.join(self._r2_telomeric_reads_clusters_directory, reads_file)), "forward")
                if self._read_type == "seperated":
                    self._r1_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._r1_reads_filepath, os.path.abspath(os.path.join(self._r2_telomeric_reads_clusters_directory, reads_file)), "forward")
        print("Done.\n")

        print("R1 subtelomeric reads saved to '" + self._r1_subtelomeric_reads_filepath + "'")
        print("R2 subtelomeric reads saved to '" + self._r2_subtelomeric_reads_filepath + "'")
        print("\n")
    
    def cluster_subtelomeric_reads(self):
        subtelomeric_reads_cluster_info_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'cluster_info', 'subtelomeric_reads'))
        os.makedirs(subtelomeric_reads_cluster_info_directory, exist_ok=True)
        
        self._r1_subtelomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(subtelomeric_reads_cluster_info_directory, 'r1_subtelomeric_reads_clusters_info.txt'))
        self._r2_subtelomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(subtelomeric_reads_cluster_info_directory, 'r2_subtelomeric_reads_clusters_info.txt'))

        # cluster forward subtelomeric reads and write results to file 
        print("Clustering subtelomeric reads in '" + self._r1_subtelomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r1_subtelomeric_reads_clusters_info_filepath + " -c " + self._r1_subtelomeric_reads_filepath + " >/dev/null 2>&1")
        self._r1_subtelomeric_reads_clusters_directory = self._generate_clustered_reads(self._r1_subtelomeric_reads_filepath, self._r1_subtelomeric_reads_clusters_info_filepath, "subtelomeric", "forward")
        print("Done.\n")

        # cluster reverse subtelomeric reads and write results to file 
        print("Clustering telomeric reads in '" + self._r2_subtelomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r2_subtelomeric_reads_clusters_info_filepath + " -c " + self._r2_subtelomeric_reads_filepath + " >/dev/null 2>&1")
        self._r2_subtelomeric_reads_clusters_directory = self._generate_clustered_reads(self._r2_subtelomeric_reads_filepath, self._r2_subtelomeric_reads_clusters_info_filepath, "subtelomeric", "reverse")
        print("Done.\n")
        
        print("R1 subtelomeric reads clustering results saved to '" + self._r1_telomeric_reads_clusters_directory)
        print("R2 subtelomeric reads clustering results saved to '" + self._r2_telomeric_reads_clusters_directory)
        print()
    
    def align_clustered_reads(self):
        if platform.system() == "Linux":
            os.system("chmod u=rwx,g=rwx,o=rwx ../lib/MUSCLE/muscle3.8.31_i86linux32")
            muscle = r"../lib/MUSCLE/muscle3.8.31_i86linux32"
        if platform.system() == "Darwin":
            os.system("chmod u=rwx,g=rwx,o=rwx ../lib/MUSCLE/muscle3.8.31_i86darwin64")
            muscle = r"../lib/MUSCLE/muscle3.8.31_i86darwin64"

        print("Aligning clustered R1 telomeric reads in '" + self._r1_telomeric_reads_clusters_directory + "'...")
        self._r1_telomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r1_telomeric_reads'))
        os.makedirs(self._r1_telomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r1_telomeric_reads_clusters_directory, self._r1_telomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Aligning clustered R2 telomeric reads in '" + self._r2_telomeric_reads_clusters_directory + "'...")
        self._r2_telomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r2_telomeric_reads'))
        os.makedirs(self._r2_telomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r2_telomeric_reads_clusters_directory, self._r2_telomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Aligning clustered R1 subtelomeric reads in '" + self._r1_subtelomeric_reads_clusters_directory + "'...")
        self._r1_subtelomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r1_subtelomeric_reads'))
        os.makedirs(self._r1_subtelomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r1_subtelomeric_reads_clusters_directory, self._r1_subtelomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Aligning clustered R2 subtelomeric reads in '" + self._r2_subtelomeric_reads_clusters_directory + "'...")
        self._r2_subtelomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r2_subtelomeric_reads'))
        os.makedirs(self._r2_subtelomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r2_subtelomeric_reads_clusters_directory, self._r2_subtelomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Alignment results for clustered R1 telomeric reads saved to '" + self._r1_telomeric_reads_clusters_alignments_directory)
        print("Alignment results for clustered R2 telomeric reads saved to '" + self._r2_telomeric_reads_clusters_alignments_directory)
        print("Alignment results for clustered R1 subtelomeric reads saved to '" + self._r1_subtelomeric_reads_clusters_alignments_directory)
        print("Alignment results for clustered R2 subtelomeric reads saved to '" + self._r2_subtelomeric_reads_clusters_alignments_directory)
        print()

    def consensus(self):
        consensus_reads_directory = os.path.abspath(os.path.join(self._out_directory, 'consensus'))
        os.makedirs(consensus_reads_directory, exist_ok=True)

        print("Generating consensus sequences for R1 telomeric reads clusters alignments...")
        self._r1_telomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r1_telomeric_reads_clusters_consensus.fasta'))
        r1_telomeric_reads_clusters_consensus_output = open(self._r1_telomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r1_telomeric_reads_clusters_alignments_directory, r1_telomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Generating consensus sequences for R2 telomeric reads clusters alignments...")
        self._r2_telomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r2_telomeric_reads_clusters_consensus.fasta'))
        r2_telomeric_reads_clusters_consensus_output = open(self._r2_telomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r2_telomeric_reads_clusters_alignments_directory, r2_telomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Generating consensus sequences for R1 subtelomeric reads clusters alignments...")
        self._r1_subtelomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r1_subtelomeric_reads_clusters_consensus.fasta'))
        r1_subtelomeric_reads_clusters_consensus_output = open(self._r1_subtelomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r1_subtelomeric_reads_clusters_alignments_directory, r1_subtelomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Generating consensus sequences for R2 subtelomeric reads clusters alignments...")
        self._r2_subtelomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r2_subtelomeric_reads_clusters_consensus.fasta'))
        r2_subtelomeric_reads_clusters_consensus_output = open(self._r2_subtelomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r2_subtelomeric_reads_clusters_alignments_directory, r2_subtelomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Consensus sequences for clustered R1 telomeric reads saved to '" + self._r1_telomeric_reads_clusters_consensus_filepath)
        print("Consensus sequences for clustered R2 telomeric reads saved to '" + self._r2_telomeric_reads_clusters_consensus_filepath)
        print("Consensus sequences for clustered R1 subtelomeric reads saved to '" + self._r1_subtelomeric_reads_clusters_consensus_filepath)
        print("Consensus sequences for clustered R2 subtelomeric reads saved to '" + self._r2_subtelomeric_reads_clusters_consensus_filepath)
        print()

    def _generate_regular_expression(self, telomeric_repeat, end):
        body = telomeric_repeat * 3
        regex = "(" + body + ")+"
        if end == "5'":
            regex = "^" + regex
        if end == "3'":
            regex = regex + "$"
        prefix = ""
        for base in reversed(telomeric_repeat[1:]):
            prefix = base + prefix
            if end == "5'":
                regex = regex + "|^" + prefix + "(" + body + ")+"
            if end == "3'":
                regex = regex + "|" + prefix + "(" + body + ")+$"
        return regex

    def _process_reads(self, input_filepath, r1_output, r2_output):
        line_counter = 0
        header_line = ""
        with open(input_filepath) as sequence_reads:
            for line in sequence_reads:
                line_counter += 1
                if line_counter % FQ_LINES_PER_READ == HEADER:
                    header_line = line
                    # replace for output to FASTA format 
                    header_line = header_line.replace("@", ">")
                if line_counter % FQ_LINES_PER_READ == SEQUENCE:
                    start_tel_exists = re.search(self._5prime_telomeric_repeat_regex, line)
                    end_tel_exists = re.search(self._3prime_telomeric_repeat_regex, line)
                    if bool(start_tel_exists) ^ bool(end_tel_exists):
                        if int(header_line.split(" ")[1].split(":")[0]) == FORWARD:
                            r1_output.write(header_line)
                            r1_output.write(line)
                        if int(header_line.split(" ")[1].split(":")[0]) == REVERSE:
                            r2_output.write(header_line)
                            r2_output.write(line)
        
    # retrieve telomeric reads based on clustering by wcdest
    def _generate_clustered_reads(self, reads_path, reads_clusters_info_path, read_type, read_direction):
        clusters = list()
        rows = list()
        row_num = 0 # keep track of row number of each cluster
        with open(reads_clusters_info_path) as clusters_file:
            # for each line in wcdest cluster file, convert to list of ints and 
            # remove newline and ending "." Adjust indexing.
            for line in clusters_file:
                row_num += 1
                line = [int(i) for i in line.strip(".\n").split(" ")]
                line = [i * FA_LINES_PER_READ for i in line]
                # filter out singleton clusters
                if len(line) > 1:
                    clusters.append(line)
                    rows.append(row_num)
        row_idx = 0
        # get all R1 or R2 telomeric reads 
        reads_file = open(reads_path, "r")
        reads_lines = reads_file.readlines()
        for cluster in clusters:
            # specify where to save FASTA files based on read direction and the 
            # type reads that are being clustered 
            if read_type == "telomeric" and read_direction == "forward":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r1_telomeric_reads'))
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r1_telomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            if read_type == "telomeric" and read_direction == "reverse":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r2_telomeric_reads'))
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r2_telomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            if read_type == "subtelomeric" and read_direction == "forward":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r1_subtelomeric_reads')) 
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r1_subtelomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            if read_type == "subtelomeric" and read_direction == "reverse":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r2_subtelomeric_reads')) 
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r2_subtelomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            # given the list of clusters with singletons filtered out, append
            # cluster reads to file 
            for i in range(len(cluster)):
                index = cluster[i]
                header_line = reads_lines[index]
                sequence_line = reads_lines[index + 1]
                # edit header to match the read direction of fastq file 
                if read_direction == "reverse":
                    header_line = header_line.replace(" 1:N", " 2:N")
                if read_direction == "forward":
                    header_line = header_line.replace(" 2:N", " 1:N")
                clustered_reads_output.write(header_line)
                clustered_reads_output.write(sequence_line)
            row_idx += 1
        return clustered_reads_output_directory

    # retrieve paired ends of a set of reads for a file with telomeric reads from
    # a corresponding FASTQ file containing the paired end
    def _retrieve_subtelomeric_reads(self, subtelomeric_reads_source_filepath, telomeric_reads_cluster_path, subtelomeric_reads_direction):
        subtelomeric_reads_headers = set()
        line_counter = 0
        # store all headers from telomeric read file
        with open(telomeric_reads_cluster_path) as telomeric_reads_cluster:
            for line in telomeric_reads_cluster:
                line_counter += 1
                if line_counter % FA_LINES_PER_READ == HEADER:
                    # edit header to match the read direction of fastq file 
                    if subtelomeric_reads_direction == "forward":
                        line = line.replace(" 2:N", " 1:N")
                    if subtelomeric_reads_direction == "reverse":
                        line = line.replace(" 1:N", " 2:N")
                    # edit header to match FASTQ format for searching original 
                    subtelomeric_reads_headers.add(line.replace(">", "@"))

        # specify where to save FASTA files based on the read type
        if subtelomeric_reads_direction == "forward":
            subtelomeric_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'subtelomeric_reads'))
            os.makedirs(subtelomeric_reads_output_directory, exist_ok=True)
            r1_subtelomeric_reads_output_filepath = os.path.abspath(os.path.join(subtelomeric_reads_output_directory, 'r1_subtelomeric_reads.fasta'))
            r1_subtelomeric_reads_output = open(r1_subtelomeric_reads_output_filepath , "a+")
        if subtelomeric_reads_direction == "reverse":
            subtelomeric_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'subtelomeric_reads'))
            os.makedirs(subtelomeric_reads_output_directory, exist_ok=True)
            r2_subtelomeric_reads_output_filepath = os.path.abspath(os.path.join(subtelomeric_reads_output_directory, 'r2_subtelomeric_reads.fasta'))
            r2_subtelomeric_reads_output = open(r2_subtelomeric_reads_output_filepath, "a+")
        line_counter = 1
        header_line = ""
        # write all paired ends to file 
        with open(subtelomeric_reads_source_filepath) as subtelomeric_reads_source:
            for line in subtelomeric_reads_source:
                if line_counter % FQ_LINES_PER_READ == HEADER:
                    header_line = line
                if line_counter % FQ_LINES_PER_READ == SEQUENCE:
                    # write sequence from FASTQ to file if its header is specified 
                    # in the list tel_headers constructed earlier 
                    if header_line in subtelomeric_reads_headers:
                        # edit header to match FASTA format
                        header_line = header_line.replace("@", ">")
                        if int(header_line.split(" ")[1].split(":")[0]) == FORWARD:
                            r1_subtelomeric_reads_output.write(header_line)
                            r1_subtelomeric_reads_output.write(line)
                        if int(header_line.split(" ")[1].split(":")[0]) == REVERSE:
                            r2_subtelomeric_reads_output.write(header_line)
                            r2_subtelomeric_reads_output.write(line)
                line_counter += 1
        if subtelomeric_reads_direction == "forward":
            return r1_subtelomeric_reads_output_filepath
        if subtelomeric_reads_direction == "reverse":
            return r2_subtelomeric_reads_output_filepath

    def _generate_alignments(self, muscle, clustered_reads_directory, alignment_output_directory):
        for cluster_reads_file in os.listdir(clustered_reads_directory):
            if cluster_reads_file.endswith(".fasta"):
                cluster_reads_filepath = os.path.abspath(os.path.join(clustered_reads_directory, cluster_reads_file))
                alignment_output_filename = cluster_reads_file#.split(".")[0] + ".aln"
                alignment_output_filepath = os.path.abspath(os.path.join(alignment_output_directory, alignment_output_filename))
                muscle_command = MuscleCommandline(muscle, input=cluster_reads_filepath, out=alignment_output_filepath)
                muscle_command()
    
    def _generate_consensus_sequences(self, clustered_reads_directory, clustered_reads_consensus_output):
        for alignment_file in os.listdir(clustered_reads_directory):
            if alignment_file.endswith(".fasta"):
                alignment_filepath = os.path.abspath(os.path.join(clustered_reads_directory, alignment_file))
                header_line = ">" + "contig" + alignment_file.split(".")[0].split("_")[-1] + "\n"
                alignment = AlignIO.read(alignment_filepath, "fasta")
                alignment_summary = AlignInfo.SummaryInfo(alignment)
                consensus_sequence = str(alignment_summary.gap_consensus(ambiguous='N')) + "\n"
                clustered_reads_consensus_output.write(header_line)
                clustered_reads_consensus_output.write(consensus_sequence)
#!/usr/bin/python

import os
import re
import platform

from Bio import AlignIO
from Bio.Seq import Seq
from Bio.Align import AlignInfo
from Bio.Align.Applications import MuscleCommandline

# define constants
HEADER = 1
SEQUENCE = 2

FORWARD = 1
REVERSE = 2

FA_LINES_PER_READ = 2
FQ_LINES_PER_READ = 4
class SequenceData:
    def __init__(self, read_type, telomeric_repeat, interleaved_read_filepath=None, r1_read_filepath=None, r2_read_filepath=None, out_directory=None):
        self._read_type = read_type

        self._interleaved_reads_filepath = interleaved_read_filepath 
        self._r1_reads_filepath = r1_read_filepath 
        self._r2_reads_filepath = r2_read_filepath

        self._out_directory = out_directory
        
        if self._out_directory == None:
            self._out_directory = os.path.abspath(os.path.join(os.path.expanduser('~'), 'terminus_output'))
            os.makedirs(os.path.dirname(self._out_directory), exist_ok=True)

        self._3prime_telomeric_repeat_regex = self._generate_regular_expression(telomeric_repeat, "3'")
        self._5prime_telomeric_repeat_regex = self._generate_regular_expression(str(Seq(telomeric_repeat).reverse_complement()), "5'")

        self._r1_telomeric_reads_filepath = None
        self._r2_telomeric_reads_filepath = None

        self._r1_telomeric_reads_clusters_info_filepath = None
        self._r2_telomeric_reads_clusters_info_filepath = None

        self._r1_telomeric_reads_clusters_directory = None
        self._r2_telomeric_reads_clusters_directory = None
        
        self._r1_telomeric_reads_clusters_alignments_directory = None
        self._r2_telomeric_reads_clusters_alignments_directory = None

        self._r1_telomeric_reads_clusters_consensus_filepath = None
        self._r2_telomeric_reads_clusters_consensus_filepath = None

        self._r1_subtelomeric_reads_filepath = None
        self._r2_subtelomeric_reads_filepath = None

        self._r1_subtelomeric_reads_clusters_info_filepath = None
        self._r2_subtelomeric_reads_clusters_info_filepath = None

        self._r1_subtelomeric_reads_clusters_directory = None
        self._r2_subtelomeric_reads_clusters_directory = None

        self._r1_subtelomeric_reads_clusters_alignments_directory = None
        self._r2_subtelomeric_reads_clusters_alignments_directory = None

        self._r1_subtelomeric_reads_clusters_consensus_filepath = None
        self._r2_subtelomeric_reads_clusters_consensus_filepath = None
    
    def get_read_type(self):
        return self._read_type

    def get_interleaved_reads_filepath(self):
        return self._interleaved_reads_filepath

    def get_r1_reads_filepath(self):
        return self._r1_reads_filepath

    def get_r2_reads_filepath(self):
        return self._r2_reads_filepath
    
    def get_out_directory(self):
        return self._out_directory
    
    def get_3prime_telomeric_repeat_regex(self):
        return self._3prime_telomeric_repeat_regex

    def get_5prime_telomeric_repeat_regex(self):
        return self._5prime_telomeric_repeat_regex
    
    def get_r1_telomeric_reads_filepath(self):
        return self._r1_telomeric_reads_filepath
    
    def get_r2_telomeric_reads_filepath(self):
        return self._r2_telomeric_reads_filepath
    
    def get_r1_telomeric_reads_clusters_info_filepath(self):
        return self._r1_telomeric_reads_clusters_info_filepath
    
    def get_r2_telomeric_reads_clusters_info_filepath(self):
        return self._r2_telomeric_reads_clusters_info_filepath

    def get_r1_telomeric_reads_clusters_directory(self):
        return self._r1_telomeric_reads_clusters_directory
    
    def get_r2_telomeric_reads_clusters_directory(self):
        return self._r2_telomeric_reads_clusters_directory

    def get_r1_telomeric_reads_clusters_alignments_directory(self):    
        return self._r1_telomeric_reads_clusters_alignments_directory
        
    def get_r2_telomeric_reads_clusters_alignments_directory(self):    
        return self._r2_telomeric_reads_clusters_alignments_directory

    def get_r1_telomeric_reads_clusters_consensus_filepath(self): 
        return self._r1_telomeric_reads_clusters_consensus_filepath
        
    def get_r2_telomeric_reads_clusters_consensus_filepath(self):
        return self._r2_telomeric_reads_clusters_consensus_filepath
    
    def get_r1_subtelomeric_reads_filepath(self):
        return self._r1_subtelomeric_reads_filepath
    
    def get_r2_subtelomeric_reads_filepath(self):
        return self._r2_subtelomeric_reads_filepath
    
    def get_r1_subtelomeric_reads_clusters_info_filepath(self):
        return self._r1_subtelomeric_reads_clusters_info_filepath
    
    def get_r2_subtelomeric_reads_clusters_info_filepath(self):
        return self._r2_subtelomeric_reads_clusters_info_filepath

    def get_r1_subtelomeric_reads_clusters_directory(self):
        return self._r1_subtelomeric_reads_clusters_directory
    
    def get_r2_subtelomeric_reads_clusters_directory(self):
        return self._r2_subtelomeric_reads_clusters_directory

    def get_r1_subtelomeric_reads_clusters_alignments_directory(self):    
        return self._r1_subtelomeric_reads_clusters_alignments_directory
        
    def get_r2_subtelomeric_reads_clusters_alignments_directory(self):    
        return self._r2_subtelomeric_reads_clusters_alignments_directory

    def get_r1_subtelomeric_reads_clusters_consensus_filepath(self): 
        return self._r1_subtelomeric_reads_clusters_consensus_filepath
        
    def get_r2_subtelomeric_reads_clusters_consensus_filepath(self):
        return self._r2_subtelomeric_reads_clusters_consensus_filepath
    
    def extract_telomeric_reads(self):
        telomeric_reads_directory = os.path.abspath(os.path.join(self._out_directory, 'telomeric_reads'))
        os.makedirs(telomeric_reads_directory, exist_ok=True)

        self._r1_telomeric_reads_filepath = os.path.abspath(os.path.join(telomeric_reads_directory, 'r1_telomeric_reads.fasta'))
        r1_telomeric_reads_output = open(self._r1_telomeric_reads_filepath, "w+")

        self._r2_telomeric_reads_filepath = os.path.abspath(os.path.join(telomeric_reads_directory, 'r2_telomeric_reads.fasta'))
        r2_telomeric_reads_output = open(self._r2_telomeric_reads_filepath, "w+")

        if self._read_type == "interleaved":
            print("Extracting telomeric reads from an interleaved FASTQ dataset...")
            self._process_reads(self._interleaved_reads_filepath, r1_telomeric_reads_output, r2_telomeric_reads_output)        
            print("Done.\n")

        if self._read_type == "seperated":
            print("Extracting telomeric reads from R1 FASTQ dataset...")
            self._process_reads(input_filepath=self._r1_reads_filepath, r1_output=r1_telomeric_reads_output, r2_output=None)
            print("Done.\n")

            print("Extracting telomeric reads from R2 FASTQ dataset...")
            self._process_reads(input_filepath=self._r2_reads_filepath, r1_output=None, r2_output=r2_telomeric_reads_output)
            print("Done.\n")       

        print("R1 telomeric reads saved to " + self._r1_telomeric_reads_filepath)
        print("R2 telomeric reads saved to " + self._r2_telomeric_reads_filepath)
        print()
 
    def cluster_telomeric_reads(self):   
        # install and make wcd for clustering
        print("Installing and preparing wcdest for clustering...")
        os.chdir("..")
        os.system("git submodule update --init --recursive >/dev/null 2>&1")
        os.chdir("./lib/wcdest/code/")
        os.system("./configure >/dev/null 2>&1")
        # NOTE: texinfo/makeinfo must be installed for wcdest to compile correctly!! 
        os.system("make clean -C ./ >/dev/null 2>&1")
        os.system("make -C ./ >/dev/null 2>&1")
        os.system("make install -C ./ >/dev/null 2>&1")
        os.chdir("../../../src")
        print("Done.\n")

        telomeric_reads_cluster_info_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'cluster_info', 'telomeric_reads'))
        os.makedirs(telomeric_reads_cluster_info_directory, exist_ok=True)
        
        self._r1_telomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(telomeric_reads_cluster_info_directory, 'r1_telomeric_reads_clusters_info.txt'))
        self._r2_telomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(telomeric_reads_cluster_info_directory, 'r2_telomeric_reads_clusters_info.txt'))

        # cluster forward telomeric reads and write results to file 
        print("Clustering telomeric reads in '" + self._r1_telomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r1_telomeric_reads_clusters_info_filepath + " -c " + self._r1_telomeric_reads_filepath + " >/dev/null 2>&1")
        self._r1_telomeric_reads_clusters_directory = self._generate_clustered_reads(self._r1_telomeric_reads_filepath, self._r1_telomeric_reads_clusters_info_filepath, "telomeric", "forward")
        print("Done.\n")

        # cluster reverse telomeric reads and write results to file 
        print("Clustering telomeric reads in '" + self._r2_telomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r2_telomeric_reads_clusters_info_filepath + " -c " + self._r2_telomeric_reads_filepath + " >/dev/null 2>&1")
        self._r2_telomeric_reads_clusters_directory = self._generate_clustered_reads(self._r2_telomeric_reads_filepath, self._r2_telomeric_reads_clusters_info_filepath, "telomeric", "reverse")
        print("Done.\n")
        
        print("R1 telomeric reads clustering results saved to '" + self._r1_telomeric_reads_clusters_directory)
        print("R2 telomeric reads clustering results saved to '" + self._r2_telomeric_reads_clusters_directory)
        print()

    def extract_subtelomeric_reads(self):
        # given forward telomeric clusters, retrieve the subtelomeric sequence reads
        #  from the reverse read file  
        print("Extracting subtelomeric reads from R2 sequence reads file for clustered R1 telomeric reads...")
        for reads_file in os.listdir(self._r1_telomeric_reads_clusters_directory):
            if reads_file.endswith(".fasta"):
                if self._read_type == "interleaved":
                    self._r2_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._interleaved_reads_filepath, os.path.abspath(os.path.join(self._r1_telomeric_reads_clusters_directory, reads_file)), "reverse")
                if self._read_type == "seperated":
                    self._r2_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._r2_reads_filepath, os.path.abspath(os.path.join(self._r1_telomeric_reads_clusters_directory, reads_file)), "reverse")
        print("Done.\n")

        # given reverse telomeric clusters, retrieve the subtelomeric sequence reads
        #  from the forward read file  
        print("Extracting subtelomeric reads from R1 sequence reads file for clustered R2 telomeric reads...")
        for reads_file in os.listdir(self._r2_telomeric_reads_clusters_directory):
            if reads_file.endswith(".fasta"): 
                if self._read_type == "interleaved":
                    self._r1_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._interleaved_reads_filepath, os.path.abspath(os.path.join(self._r2_telomeric_reads_clusters_directory, reads_file)), "forward")
                if self._read_type == "seperated":
                    self._r1_subtelomeric_reads_filepath = self._retrieve_subtelomeric_reads(self._r1_reads_filepath, os.path.abspath(os.path.join(self._r2_telomeric_reads_clusters_directory, reads_file)), "forward")
        print("Done.\n")

        print("R1 subtelomeric reads saved to '" + self._r1_subtelomeric_reads_filepath + "'")
        print("R2 subtelomeric reads saved to '" + self._r2_subtelomeric_reads_filepath + "'")
        print("\n")
    
    def cluster_subtelomeric_reads(self):
        subtelomeric_reads_cluster_info_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'cluster_info', 'subtelomeric_reads'))
        os.makedirs(subtelomeric_reads_cluster_info_directory, exist_ok=True)
        
        self._r1_subtelomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(subtelomeric_reads_cluster_info_directory, 'r1_subtelomeric_reads_clusters_info.txt'))
        self._r2_subtelomeric_reads_clusters_info_filepath = os.path.abspath(os.path.join(subtelomeric_reads_cluster_info_directory, 'r2_subtelomeric_reads_clusters_info.txt'))

        # cluster forward subtelomeric reads and write results to file 
        print("Clustering subtelomeric reads in '" + self._r1_subtelomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r1_subtelomeric_reads_clusters_info_filepath + " -c " + self._r1_subtelomeric_reads_filepath + " >/dev/null 2>&1")
        self._r1_subtelomeric_reads_clusters_directory = self._generate_clustered_reads(self._r1_subtelomeric_reads_filepath, self._r1_subtelomeric_reads_clusters_info_filepath, "subtelomeric", "forward")
        print("Done.\n")

        # cluster reverse subtelomeric reads and write results to file 
        print("Clustering telomeric reads in '" + self._r2_subtelomeric_reads_filepath + "'...")
        os.system("../lib/wcdest/code/src/wcd -o " + self._r2_subtelomeric_reads_clusters_info_filepath + " -c " + self._r2_subtelomeric_reads_filepath + " >/dev/null 2>&1")
        self._r2_subtelomeric_reads_clusters_directory = self._generate_clustered_reads(self._r2_subtelomeric_reads_filepath, self._r2_subtelomeric_reads_clusters_info_filepath, "subtelomeric", "reverse")
        print("Done.\n")
        
        print("R1 subtelomeric reads clustering results saved to '" + self._r1_telomeric_reads_clusters_directory)
        print("R2 subtelomeric reads clustering results saved to '" + self._r2_telomeric_reads_clusters_directory)
        print()
    
    def align_clustered_reads(self):
        if platform.system() == "Linux":
            os.system("chmod u=rwx,g=rwx,o=rwx ../lib/MUSCLE/muscle3.8.31_i86linux32")
            muscle = r"../lib/MUSCLE/muscle3.8.31_i86linux32"
        if platform.system() == "Darwin":
            os.system("chmod u=rwx,g=rwx,o=rwx ../lib/MUSCLE/muscle3.8.31_i86darwin64")
            muscle = r"../lib/MUSCLE/muscle3.8.31_i86darwin64"

        print("Aligning clustered R1 telomeric reads in '" + self._r1_telomeric_reads_clusters_directory + "'...")
        self._r1_telomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r1_telomeric_reads'))
        os.makedirs(self._r1_telomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r1_telomeric_reads_clusters_directory, self._r1_telomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Aligning clustered R2 telomeric reads in '" + self._r2_telomeric_reads_clusters_directory + "'...")
        self._r2_telomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r2_telomeric_reads'))
        os.makedirs(self._r2_telomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r2_telomeric_reads_clusters_directory, self._r2_telomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Aligning clustered R1 subtelomeric reads in '" + self._r1_subtelomeric_reads_clusters_directory + "'...")
        self._r1_subtelomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r1_subtelomeric_reads'))
        os.makedirs(self._r1_subtelomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r1_subtelomeric_reads_clusters_directory, self._r1_subtelomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Aligning clustered R2 subtelomeric reads in '" + self._r2_subtelomeric_reads_clusters_directory + "'...")
        self._r2_subtelomeric_reads_clusters_alignments_directory = os.path.abspath(os.path.join(self._out_directory, 'alignments', 'r2_subtelomeric_reads'))
        os.makedirs(self._r2_subtelomeric_reads_clusters_alignments_directory, exist_ok=True)
        self._generate_alignments(muscle, self._r2_subtelomeric_reads_clusters_directory, self._r2_subtelomeric_reads_clusters_alignments_directory)
        print("Done.\n")

        print("Alignment results for clustered R1 telomeric reads saved to '" + self._r1_telomeric_reads_clusters_alignments_directory)
        print("Alignment results for clustered R2 telomeric reads saved to '" + self._r2_telomeric_reads_clusters_alignments_directory)
        print("Alignment results for clustered R1 subtelomeric reads saved to '" + self._r1_subtelomeric_reads_clusters_alignments_directory)
        print("Alignment results for clustered R2 subtelomeric reads saved to '" + self._r2_subtelomeric_reads_clusters_alignments_directory)
        print()

    def consensus(self):
        consensus_reads_directory = os.path.abspath(os.path.join(self._out_directory, 'consensus'))
        os.makedirs(consensus_reads_directory, exist_ok=True)

        print("Generating consensus sequences for R1 telomeric reads clusters alignments...")
        self._r1_telomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r1_telomeric_reads_clusters_consensus.fasta'))
        r1_telomeric_reads_clusters_consensus_output = open(self._r1_telomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r1_telomeric_reads_clusters_alignments_directory, r1_telomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Generating consensus sequences for R2 telomeric reads clusters alignments...")
        self._r2_telomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r2_telomeric_reads_clusters_consensus.fasta'))
        r2_telomeric_reads_clusters_consensus_output = open(self._r2_telomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r2_telomeric_reads_clusters_alignments_directory, r2_telomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Generating consensus sequences for R1 subtelomeric reads clusters alignments...")
        self._r1_subtelomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r1_subtelomeric_reads_clusters_consensus.fasta'))
        r1_subtelomeric_reads_clusters_consensus_output = open(self._r1_subtelomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r1_subtelomeric_reads_clusters_alignments_directory, r1_subtelomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Generating consensus sequences for R2 subtelomeric reads clusters alignments...")
        self._r2_subtelomeric_reads_clusters_consensus_filepath = os.path.abspath(os.path.join(consensus_reads_directory, 'r2_subtelomeric_reads_clusters_consensus.fasta'))
        r2_subtelomeric_reads_clusters_consensus_output = open(self._r2_subtelomeric_reads_clusters_consensus_filepath, "w+")
        self._generate_consensus_sequences(self._r2_subtelomeric_reads_clusters_alignments_directory, r2_subtelomeric_reads_clusters_consensus_output)
        print("Done.\n")

        print("Consensus sequences for clustered R1 telomeric reads saved to '" + self._r1_telomeric_reads_clusters_consensus_filepath)
        print("Consensus sequences for clustered R2 telomeric reads saved to '" + self._r2_telomeric_reads_clusters_consensus_filepath)
        print("Consensus sequences for clustered R1 subtelomeric reads saved to '" + self._r1_subtelomeric_reads_clusters_consensus_filepath)
        print("Consensus sequences for clustered R2 subtelomeric reads saved to '" + self._r2_subtelomeric_reads_clusters_consensus_filepath)
        print()

    def _generate_regular_expression(self, telomeric_repeat, end):
        body = telomeric_repeat * 3
        regex = "(" + body + ")+"
        if end == "5'":
            regex = "^" + regex
        if end == "3'":
            regex = regex + "$"
        prefix = ""
        for base in reversed(telomeric_repeat[1:]):
            prefix = base + prefix
            if end == "5'":
                regex = regex + "|^" + prefix + "(" + body + ")+"
            if end == "3'":
                regex = regex + "|" + prefix + "(" + body + ")+$"
        return regex

    def _process_reads(self, input_filepath, r1_output, r2_output):
        line_counter = 0
        header_line = ""
        with open(input_filepath) as sequence_reads:
            for line in sequence_reads:
                line_counter += 1
                if line_counter % FQ_LINES_PER_READ == HEADER:
                    header_line = line
                    # replace for output to FASTA format 
                    header_line = header_line.replace("@", ">")
                if line_counter % FQ_LINES_PER_READ == SEQUENCE:
                    start_tel_exists = re.search(self._5prime_telomeric_repeat_regex, line)
                    end_tel_exists = re.search(self._3prime_telomeric_repeat_regex, line)
                    if bool(start_tel_exists) ^ bool(end_tel_exists):
                        if int(header_line.split(" ")[1].split(":")[0]) == FORWARD:
                            r1_output.write(header_line)
                            r1_output.write(line)
                        if int(header_line.split(" ")[1].split(":")[0]) == REVERSE:
                            r2_output.write(header_line)
                            r2_output.write(line)
        
    # retrieve telomeric reads based on clustering by wcdest
    def _generate_clustered_reads(self, reads_path, reads_clusters_info_path, read_type, read_direction):
        clusters = list()
        rows = list()
        row_num = 0 # keep track of row number of each cluster
        with open(reads_clusters_info_path) as clusters_file:
            # for each line in wcdest cluster file, convert to list of ints and 
            # remove newline and ending "." Adjust indexing.
            for line in clusters_file:
                row_num += 1
                line = [int(i) for i in line.strip(".\n").split(" ")]
                line = [i * FA_LINES_PER_READ for i in line]
                # filter out singleton clusters
                if len(line) > 1:
                    clusters.append(line)
                    rows.append(row_num)
        row_idx = 0
        # get all R1 or R2 telomeric reads 
        reads_file = open(reads_path, "r")
        reads_lines = reads_file.readlines()
        for cluster in clusters:
            # specify where to save FASTA files based on read direction and the 
            # type reads that are being clustered 
            if read_type == "telomeric" and read_direction == "forward":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r1_telomeric_reads'))
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r1_telomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            if read_type == "telomeric" and read_direction == "reverse":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r2_telomeric_reads'))
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r2_telomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            if read_type == "subtelomeric" and read_direction == "forward":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r1_subtelomeric_reads')) 
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r1_subtelomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            if read_type == "subtelomeric" and read_direction == "reverse":
                clustered_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'clusters', 'r2_subtelomeric_reads')) 
                os.makedirs(clustered_reads_output_directory, exist_ok=True)
                clustered_reads_output_filepath = os.path.abspath(os.path.join(clustered_reads_output_directory, 'r2_subtelomeric_reads_cluster_' + str(rows[row_idx]) + '.fasta'))
                clustered_reads_output = open(clustered_reads_output_filepath, "w+")
            # given the list of clusters with singletons filtered out, append
            # cluster reads to file 
            for i in range(len(cluster)):
                index = cluster[i]
                header_line = reads_lines[index]
                sequence_line = reads_lines[index + 1]
                # edit header to match the read direction of fastq file 
                if read_direction == "reverse":
                    header_line = header_line.replace(" 1:N", " 2:N")
                if read_direction == "forward":
                    header_line = header_line.replace(" 2:N", " 1:N")
                clustered_reads_output.write(header_line)
                clustered_reads_output.write(sequence_line)
            row_idx += 1
        return clustered_reads_output_directory

    # retrieve paired ends of a set of reads for a file with telomeric reads from
    # a corresponding FASTQ file containing the paired end
    def _retrieve_subtelomeric_reads(self, subtelomeric_reads_source_filepath, telomeric_reads_cluster_path, subtelomeric_reads_direction):
        subtelomeric_reads_headers = set()
        line_counter = 0
        # store all headers from telomeric read file
        with open(telomeric_reads_cluster_path) as telomeric_reads_cluster:
            for line in telomeric_reads_cluster:
                line_counter += 1
                if line_counter % FA_LINES_PER_READ == HEADER:
                    # edit header to match the read direction of fastq file 
                    if subtelomeric_reads_direction == "forward":
                        line = line.replace(" 2:N", " 1:N")
                    if subtelomeric_reads_direction == "reverse":
                        line = line.replace(" 1:N", " 2:N")
                    # edit header to match FASTQ format for searching original 
                    subtelomeric_reads_headers.add(line.replace(">", "@"))

        # specify where to save FASTA files based on the read type
        if subtelomeric_reads_direction == "forward":
            subtelomeric_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'subtelomeric_reads'))
            os.makedirs(subtelomeric_reads_output_directory, exist_ok=True)
            r1_subtelomeric_reads_output_filepath = os.path.abspath(os.path.join(subtelomeric_reads_output_directory, 'r1_subtelomeric_reads.fasta'))
            r1_subtelomeric_reads_output = open(r1_subtelomeric_reads_output_filepath , "a+")
        if subtelomeric_reads_direction == "reverse":
            subtelomeric_reads_output_directory = os.path.abspath(os.path.join(self._out_directory, 'subtelomeric_reads'))
            os.makedirs(subtelomeric_reads_output_directory, exist_ok=True)
            r2_subtelomeric_reads_output_filepath = os.path.abspath(os.path.join(subtelomeric_reads_output_directory, 'r2_subtelomeric_reads.fasta'))
            r2_subtelomeric_reads_output = open(r2_subtelomeric_reads_output_filepath, "a+")
        line_counter = 1
        header_line = ""
        # write all paired ends to file 
        with open(subtelomeric_reads_source_filepath) as subtelomeric_reads_source:
            for line in subtelomeric_reads_source:
                if line_counter % FQ_LINES_PER_READ == HEADER:
                    header_line = line
                if line_counter % FQ_LINES_PER_READ == SEQUENCE:
                    # write sequence from FASTQ to file if its header is specified 
                    # in the list tel_headers constructed earlier 
                    if header_line in subtelomeric_reads_headers:
                        # edit header to match FASTA format
                        header_line = header_line.replace("@", ">")
                        if int(header_line.split(" ")[1].split(":")[0]) == FORWARD:
                            r1_subtelomeric_reads_output.write(header_line)
                            r1_subtelomeric_reads_output.write(line)
                        if int(header_line.split(" ")[1].split(":")[0]) == REVERSE:
                            r2_subtelomeric_reads_output.write(header_line)
                            r2_subtelomeric_reads_output.write(line)
                line_counter += 1
        if subtelomeric_reads_direction == "forward":
            return r1_subtelomeric_reads_output_filepath
        if subtelomeric_reads_direction == "reverse":
            return r2_subtelomeric_reads_output_filepath

    def _generate_alignments(self, muscle, clustered_reads_directory, alignment_output_directory):
        for cluster_reads_file in os.listdir(clustered_reads_directory):
            if cluster_reads_file.endswith(".fasta"):
                cluster_reads_filepath = os.path.abspath(os.path.join(clustered_reads_directory, cluster_reads_file))
                alignment_output_filename = cluster_reads_file#.split(".")[0] + ".aln"
                alignment_output_filepath = os.path.abspath(os.path.join(alignment_output_directory, alignment_output_filename))
                muscle_command = MuscleCommandline(muscle, input=cluster_reads_filepath, out=alignment_output_filepath)
                muscle_command()
    
    def _generate_consensus_sequences(self, clustered_reads_directory, clustered_reads_consensus_output):
        for alignment_file in os.listdir(clustered_reads_directory):
            if alignment_file.endswith(".fasta"):
                alignment_filepath = os.path.abspath(os.path.join(clustered_reads_directory, alignment_file))
                header_line = ">" + "contig" + alignment_file.split(".")[0].split("_")[-1] + "\n"
                alignment = AlignIO.read(alignment_filepath, "fasta")
                alignment_summary = AlignInfo.SummaryInfo(alignment)
                consensus_sequence = str(alignment_summary.gap_consensus(ambiguous='N')) + "\n"
                clustered_reads_consensus_output.write(header_line)
                clustered_reads_consensus_output.write(consensus_sequence)
#!/usr/bin/python # define constants # install and make wcd for clustering # NOTE: texinfo/makeinfo must be installed for wcdest to compile correctly!! # cluster forward telomeric reads and write results to file # cluster reverse telomeric reads and write results to file # given forward telomeric clusters, retrieve the subtelomeric sequence reads #  from the reverse read file # given reverse telomeric clusters, retrieve the subtelomeric sequence reads #  from the forward read file # cluster forward subtelomeric reads and write results to file # cluster reverse subtelomeric reads and write results to file # replace for output to FASTA format # retrieve telomeric reads based on clustering by wcdest # keep track of row number of each cluster # for each line in wcdest cluster file, convert to list of ints and # remove newline and ending "." Adjust indexing. # filter out singleton clusters # get all R1 or R2 telomeric reads # specify where to save FASTA files based on read direction and the # type reads that are being clustered # given the list of clusters with singletons filtered out, append # cluster reads to file # edit header to match the read direction of fastq file # retrieve paired ends of a set of reads for a file with telomeric reads from # a corresponding FASTQ file containing the paired end # store all headers from telomeric read file # edit header to match the read direction of fastq file # edit header to match FASTQ format for searching original # specify where to save FASTA files based on the read type # write all paired ends to file # write sequence from FASTQ to file if its header is specified # in the list tel_headers constructed earlier # edit header to match FASTA format #.split(".")[0] + ".aln"