main.cpp

/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/* 
 * File:   main.cpp
 * Author: madiha
 *
 * Created on October 28 2020, 3:31 PM
 */

#include <cstdlib>
#include "utils.h"
#include "global.h"
#include "bySoftclip.h"
#include "mis_assembly_chimera.h"
#include "alignment.h"
#include <sstream>
#include <iostream>
#include <ctime>
#include <unistd.h>
#include <boost/filesystem.hpp>
#include "filterSamToFastq.h"
#include "byRI.h"
using namespace boost;
using namespace std;

void help() {

    cerr << endl << "###############################################################################################################################" << endl << endl;

    cerr << "Program: ROAST (Reference free Optimization of Assembled Supertranscriptomes)" << endl << "Version: 1.0.0 (using bamtools 2.4.0 and boost c++ libraries)" << endl << endl;

    cerr << "Usage: ROAST --fastq_1 fastq1_filename --fastq_2 fastq2_filename <Reference type> --threads INT [Parameters]" << endl << endl;

    cerr << "Reference type:" << endl;

    cerr << "   --supertranscript_assembly <STR> \t Provide SuperTranscript fasta for supertranscript assembly improvement" << endl;

    cerr << "   --trinity_assembly <STR>  \t\t Provide transcript fasta for raw Trinity assembly improvement" << endl;

    cerr << "   --generate_assembly \t\t\t Set and export the environmental variable TRINITY_HOME to point Trinity installation folder to generate De novo Transcriptome assembly and improve output" << endl << endl;


    cerr << "Parameters:" << endl;

    cerr << "   --output_dir <STR>  \t\t\t Path for output directories, default folder of input reference sequence" << endl;

    cerr << "   --inner_itr <INT> \t\t\t Number of inner iterations threshold, default 30" << endl;

    cerr << "   --outer_itr <INT> \t\t\t Number of outer Iterations threshold, default 100" << endl;

    cerr << "   --improvment_TH <INT> \t\t Keep improving until number of improved contigs meet threshold, default 1" << endl;

    cerr << "   --min_extended_contigs <INT> \t Keep extending until number of extended contigs meet threshold, default 1" << endl;

    cerr << "   --max_memory_TRINITY <INT> \t\t Maximum memory for TRINITY, default 20" << endl;

    cerr << "   --threads <INT> \t\t\t Number of threads, default 8" << endl;

    cerr << "   --threadsForSamSort <INT> \t\t Number of threads for samtools sort, default 2" << endl;

    cerr << "   --memForSamSort <INT> \t\t Memory for samtools sort, default 768M" << endl;

    cerr << "   --complete_cleanup \t\t\t Complete cleanup of intermediate files, default partial cleanup" << endl;

    cerr << "   --nochange_header \t\t\t Don't change header of final improved assembly file, default change header" << endl;

    cerr << "   --cdHitEST <INT> \t\t\t 0 for no CD-HIT-EST, 1 for CD-HIT-EST only in the start, 2 for CD-HIT-EST after every iteration, default 1" << endl;

    cerr << "   --mapping_quality <INT> \t\t Minimum mapping quality to filter Bam file, default 20" << endl << endl;

    cerr << "   --sc_support_TH <FLOAT> \t\t Softclips support for a position for Incomplete and fragmented contigs and false chimera process, default 0.75" << endl;

    cerr << "   --min_SCs_reads <INT> \t\t Minimum number of softclips reads for extension of Incomplete contigs, default 3" << endl;

    cerr << "   --sc_cons_len_TH <INT> \t\t Minimum length of consensus sequence generated by softclips to BLAST for fragmented contigs Identification, default 10" << endl;

    cerr << "   --discard_contig_corner_len_TH <INT>  Minimum number of bases allowed to discard from the corner of contigs while finding overlapped edges, default 10" << endl;

    cerr << "   --max_allowed_gaps <INT> \t\t Maximum gaps allowed to find overlaps between the corners of the two contigs, default 0" << endl;

    cerr << "   --edge_boundary <INT> \t\t Terminal region for mapped reads to be considered for un-mapped and distantly mapped reads, default 2x of read length" << endl;

    cerr << "   --min_unmapped_reads <INT> \t\t Minimum number of unmapped reads to generate CAP3 assembly for incomplete contigs extension, default 5" << endl;

    cerr << "   --min_distMapped_reads <INT> \t Minimum number of distantly mapped reads to consider as read island for merging of fragmented contigs, default 3" << endl;

    cerr << "   --contig_boundary   \t\t\t Check overlap between two fragmented contigs within 5% read length of the contig boundary, default read island boundary" << endl;

    cerr << "   --min_allowed_unmapped_ext <INT> \t Minimum number of bases to be considered for valid extension using unmapped reads, default 50% of read length" << endl;

    cerr << "   --one_side_allowed_SCs_RI <INT> \t Maximum % of softclips of the total read length allowed at one side of read to consider it for Read island, defaultf 25" << endl;

    cerr << "   --each_side_allowed_SCs_RI <INT> \t Maximum % of softclips of the total read length allowed at both sides of read to consider it for Read island, default 12" << endl;

    cerr << "   --sc_start_pos_from_terminus <INT> \t For terminal softclips extraction, define starting position to consider, default 25" << endl;

    cerr << "   --min_overlap_TH <INT> \t\t Minimum length of overlapped sequence to consider for merging fragmented contigs and CAP3 assemblies, default 20 bases" << endl;

    cerr << "   --win_size <INT> \t\t\t Window size to detect gradual coverage change, default 2X read_length " << endl;
    
    cerr << "   --win_diff_TH <INT> \t\t\t Maximum average coverage change threshold between two consecutive windows, default off " << endl;

    cerr << "   --coverage_drop_TH <FLOAT> \t\t Maximum coverage ratio between two consecutive positions to process for false chimera identification process, default 0.2" << endl;

    cerr << "   --st_end_boundary <INT> \t\t Consider coverage change within start and end boundary of the contig, default one and half of read length " << endl;

    cerr << "   --blast_score_TH <INT> \t\t BLAST hit identity and coverage score to find overlap between fragmented contigs and identify mis-assembly/false chimera, default 90" << endl;

    cerr << "   --ignore_short_seq <INT> \t\t Minimum length of the contig to remove from the assembly, default200" << endl;

    cerr << "   --insert_Ns <INT> \t\t\t Number of Ns to insert between contig and its CAP3 assembly in the absence of overlap, default 5" << endl << endl;

    cerr << "Terminate ROAST process:" << endl;
    
    cerr << "   To stop ROAST properly before completion of default iterations place empty file named 'stop.txt' in the folder 'intermediate_Improved_assemblies'." << endl;
  
    cerr << endl << "###############################################################################################################################" << endl << endl;
}

std::string getexepath() {
    char result[ PATH_MAX ];
    ssize_t count = readlink("/proc/self/exe", result, PATH_MAX);
    return std::string(result, (count > 0) ? count : 0);
}

bool isNumber(string s) {
    for (int i = 0; i < s.length(); i++)
        if (isdigit(s[i]) == false)
            return false;

    return true;
}

int main(int argc, char** argv) {
    if (argc < 2) {
        cout << "No arguments found" << endl;
        help();
        exit(0);
    }
    utils utils;

    //samtools --version | cut -d' ' -f2

    string samtools_version = "samtools_version.txt";
    string samtools_version_check = "samtools --version | cut -d ' ' -f2  > " + samtools_version;
    std::system(samtools_version_check.c_str());

    ifstream samtools;
    samtools.open(samtools_version.c_str());
    string line;
    stringstream v;
    int version;
    vector <string> temp;

    while (getline(samtools, line)) {
        utils.str_split(line, temp, ".");
        stringstream v(temp[1]);
        // cout << v.str() << endl;
        v >> version;

        if (version < 9) {
            cerr << "Please update your samtools to >= 1.9." << endl;
            exit(0);
        }
    }
    samtools.close();
    utils.remove_file(samtools_version);

    vector <string> args(argv, argv + argc);

    std::string s;

    for (std::vector<std::string>::const_iterator i = args.begin(); i != args.end(); ++i)
        s += " " + *i;

    cout << s << endl << endl; // print user command in the start

    string fastq_first = "";
    string fastq_sec = "";
    string fasta = "";

    string ref_ST_CD_fasta, ST_fasta, CDhitEST_command, temp_assembly;
    size_t find;

    string get_ful_path = getexepath();
    size_t exe = get_ful_path.find_last_of("/\\");
    exe_path = get_ful_path.substr(0, exe + 1);


    string ROAST_extend_contigs = exe_path + "ROAST_extendContigs";

    if (boost::filesystem::exists(ROAST_extend_contigs)) // does filePath actually exist?
    {
        //file found
    } else {
        cerr << "Please compile ROAST_extend_contigs.cpp file first" << endl;
        exit(0);
    }

    string ROAST_extendContigs_SCs = exe_path + "ROAST_extendContigs_SCs";

    if (boost::filesystem::exists(ROAST_extendContigs_SCs)) // does filePath actually exist?
    {
        //file found
    } else {
        cerr << "Please compile ROAST_extendContigs_SCs.cpp file first" << endl;
        exit(0);
    }

    string ROAST_mergeContigs_SCs = exe_path + "ROAST_mergeContigs_SCs";

    if (boost::filesystem::exists(ROAST_mergeContigs_SCs)) // does filePath actually exist?
    {
        //file found
    } else {
        cerr << "Please compile ROAST_mergeContigs_SCs.cpp file first" << endl;
        exit(0);
    }

   /* string ROAST_mergeContigs_SCs_R2 = exe_path + "ROAST_mergeContigs_SCs_R2";

    if (boost::filesystem::exists(ROAST_mergeContigs_SCs_R2)) // does filePath actually exist?
    {
        //file found
    } else {
        cerr << "Please compile ROAST_mergeContigs_SCs_R2.cpp file first" << endl;
        exit(0);
    }*/

    for (int i = 1; i < argc; i++) {

        if (args[i] == "--supertranscript_assembly") {

            ST_fasta = args[i + 1];
            //fasta = ref_ST_CD_fasta;//;.substr(0, ref_ST_CD_fasta.size() - 6) + "_0.fasta";

            if (boost::filesystem::exists(ST_fasta)) // does filePath actually exist?
            {
                //file found
            } else {
                cerr << "supertranscript_assembly doesn't have valid file path or name" << endl;
                exit(0);
            }

            fasta = ST_fasta;

        } else if (args[i] == "--trinity_assembly") {

            string trinity_assembly = args[i + 1];

            if (boost::filesystem::exists(trinity_assembly)) // does filePath actually exist?
            {
                //file found
            } else {
                cerr << "trinity_assembly doesn't have valid file path or name" << endl;
                exit(0);
            }

            //TRINITY_HOME/Analysis/OuterTranscripts/Trinity_gene_splice_modeler.py --trinity_fasta AL.fasta --out_prefix AL.ST.fasta
            find = trinity_assembly.find_last_of(".");
            ST_fasta = trinity_assembly.substr(0, find) + ".ST";
            string trinitySuperTranscript_command = "python " + exe_path + "external_tools/SuperTranscripts/Trinity_gene_splice_modeler.py --trinity_fasta " + trinity_assembly + " --out_prefix " + ST_fasta + "> /dev/null 2>&1";
           // cout << trinitySuperTranscript_command << endl;
            std::system(trinitySuperTranscript_command.c_str());

            ST_fasta = ST_fasta + ".fasta";

            if (boost::filesystem::exists(ST_fasta)) // does supertranscript file exist?
            {
                //nothing
            } else {
                cerr << endl << "Couldn't generate SuperTranscript from given Trinity assembly. Please inspect the error or provide generated ST for improvement...." << endl;
                exit(0);
            }
            fasta = ST_fasta;

        } else if (args[i] == "--fastq_1") {

            fastq_first = args[i + 1];

            if (boost::filesystem::exists(fastq_first)) // does filePath actually exist?
            {
                //file found
            } else {
                cerr << endl << "fastq_1 doesn't have valid file path or name." << endl;
                exit(0);
            }

        } else if (args[i] == "--fastq_2") {

            fastq_sec = args[i + 1];
            if (boost::filesystem::exists(fastq_sec)) // does filePath actually exist?
            {
                //file found
            } else {
                cerr << endl << "fastq_2 doesn't have valid file path or name." << endl;
                exit(0);
            }

        } else if (args[i] == "--generate_assembly") {

            generate_assembly = true;

        } else if (args[i] == "--max_memory_TRINITY") {

            if (isNumber(args[i + 1]))
                max_memory_TRINITY = args[i + 1];
            else {
                cerr << "Provided maximum memory is not valid." << endl;
                return 1;
            }

        } else if (args[i] == "--threads") {

            allowed_threads = atoi(args[i + 1].c_str());

            if (isNumber(args[i + 1]) && allowed_threads > 0)
                threads = args[i + 1];
            else {
                cerr << "Provided threads are not valid." << endl;
                return 1;
            }

        } else if (args[i] == "--threadsForSamSort") {

            threadsForSamSort = atoi(args[i + 1].c_str());

            if (isNumber(args[i + 1]) && threadsForSamSort > 0)
                threadsForSamSortS = args[i + 1];
            else {
                cerr << "Provided threads for samtools sort are not valid." << endl;
                return 1;
            }

        } else if (args[i] == "--memForSamSort") {

            threadsForSamSort = atoi(args[i + 1].c_str());

            if (isNumber(args[i + 1]) && MemForSamSort > 0)
                MemForSamSortS = args[i + 1];
            else {
                cerr << "Provided Memory for samtools sort are not valid." << endl;
                return 1;
            }

        } else if (args[i] == "--inner_itr") {

            inner_itr_TH = atoi(args[i + 1].c_str());

            if (inner_itr_TH < 0) {
                cerr << "Number of Inner iterations must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--improvment_TH") {

            improvment_TH = atoi(args[i + 1].c_str());

            if (improvment_TH < 0) {
                cerr << "Number of minimum improved contigs must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--min_extended_contigs") {

            min_extended_contigs = atoi(args[i + 1].c_str());

            if (min_extended_contigs < 0) {
                cerr << "Number of minimum extended contigs must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--outer_itr") {

            outer_itr_TH = atoi(args[i + 1].c_str());

            if (outer_itr_TH < 0) {
                cerr << "Number of Outer iterations must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--complete_cleanup") {

            complete_cleanup = true;

        } else if (args[i] == "--trnasRate_score") {

            transRate_score = true;

        } else if (args[i] == "--nochange_header") {

            change_header = false;

        } else if (args[i] == "--output_dir") {

            out_directory = args[i + 1];

            if (boost::filesystem::exists(out_directory)) // does filePath actually exist?
            {
                //file found
            } else { // make directory if not exist already
                string mkdir = "mkdir " + out_directory;
                std::system(mkdir.c_str());
                //exit(0);
            }

        } else if (args[i] == "--cdHitEST") {

              cdhitest = atoi(args[i + 1].c_str());

            if (cdhitest < 0 || cdhitest > 2) {
                cerr << "CD-HIT-EST option can be 0 to stop, 1 to run only in the start and 2 to run after super iteration" << endl;
                return 1;
            }

        } else if (args[i] == "--sc_support_TH") {

            SC_support_TH = atoi(args[i + 1].c_str());

            if (SC_support_TH > 0 && SC_support_TH < 1) {
                SC_support_TH = SC_support_TH * 100;
            } else {
                cerr << "Threshold for softclips support for a position must be between 0.0 and 1.0." << endl;
                return 1;
            }

        } else if (args[i] == "--mapping_quality") {

            mapping_quality_TH = atoi(args[i + 1].c_str());

            if (mapping_quality_TH < 0) {
                cerr << "Mapping quality to filter Bam file must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--min_SCs_reads") {

            min_sc_reads = atoi(args[i + 1].c_str());

            if (min_sc_reads < 0) {
                cerr << "Minimum number of softclips reads for extension of Incomplete contigs must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--sc_cons_len_TH") {

            SC_cons_len_TH = atoi(args[i + 1].c_str());

            if (SC_cons_len_TH < 0) {
                cerr << "Minimum length of consensus sequence generated by Softclips must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--discard_contig_corner_len_TH") {

            discard_contig_cor_len = atoi(args[i + 1].c_str());

            if (discard_contig_cor_len < 0) {
                cerr << "Minimum number of bases to discard at the corner of the contig for overlap search must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--max_allowed_gaps") {

            max_allowed_gaps = atoi(args[i + 1].c_str());

            if (max_allowed_gaps < 0) {
                cerr << "Maximum number of gaps allowed for BLAST overlap search must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--min_unmapped_reads") { // 

            min_unmapped_reads_CAP3 = atoi(args[i + 1].c_str());

            if (min_unmapped_reads_CAP3 < 0) {
                cerr << "Minimum number of unmapped reads for CAP3 assembly generation must be non-negative." << endl;
                return 1;
            }
        } else if (args[i] == "--min_distMapped_reads") {

            min_distMapped_reads = atoi(args[i + 1].c_str());

            if (min_distMapped_reads < 0) {
                cerr << "Minimum number of distantly mapped reads to consider as Read Island must be non-negative." << endl;
                return 1;
            }
        } else if (args[i] == "--min_allowed_unmapped_ext") {

            max_allowed_gaps = atoi(args[i + 1].c_str());

            if (min_CAP3_ext < 0) {
                cerr << "Minimum number of bases to be considered as valid extension using unmapped reads must be non-negative." << endl;
                return 1;
            }
        } else if (args[i] == "--contig_boundary") {

            contig_boundary = true;

        } else if (args[i] == "--one_side_allowed_SCs_RI") {

            one_side_allowed_SCs_RI = atoi(args[i + 1].c_str());

            if (one_side_allowed_SCs_RI < 0) {
                cerr << "Maximum percentage of one_side_allowed_SCs_RI must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--each_side_allowed_SCs_RI") {

            each_side_allowed_SCs_RI = atoi(args[i + 1].c_str());

            if (each_side_allowed_SCs_RI < 0) {
                cerr << "Maximum percentage of each_side_allowed_SCs_RI must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--sc_start_pos_from_terminus") {

            sc_pos_from_corner = atoi(args[i + 1].c_str());

            if (sc_pos_from_corner < 0) {
                cerr << "For terminal Softclips extraction position must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--edge_boundary") {

            left_edge_boundary = atoi(args[i + 1].c_str());

            if (left_edge_boundary < 0) {
                cerr << "Terminal region for mapped reads to be considered for un-mapped and distantly mapped  must be non-negative." << endl;
                return 1;
            }


        } else if (args[i] == "--min_overlap_TH") {

            min_overlap_TH = atoi(args[i + 1].c_str());

            if (min_overlap_TH < 0) {
                cerr << "Minimum length of overlapped sequence must be non-negative." << endl;
                return 1;
            }

        } else if (args[i] == "--st_end_boundary") {

            st_end_th = atoi(args[i + 1].c_str());

            if (st_end_th <= 0) {
                cerr << "Start and end boundary of a contig to consider coverage change must be non-negative ." << endl;
                return 1;
            }

        } else if (args[i] == "--coverage_drop_TH") {

            coverage_drop = atof(args[i + 1].c_str());      

            if (coverage_drop < 0 || coverage_drop > 1) {
                cerr << "Coverage drop threshold between two positions must be between 0.0 and 1.0." << endl;
                return 1;
            }

        } else if (args[i] == "--win_size") {

            win_size = atof(args[i + 1].c_str());

            if (win_size < 0 ) {
                cerr << "Window size of bases to identify gradual coverage change must be  anon-negative number" << endl;
                return 1;
            }

        } else if (args[i] == "--win_diff_TH") {

            win_diff_TH = atof(args[i + 1].c_str());

            if (win_diff_TH < 0 || win_diff_TH > 1) {
                cerr << "Average coverage change threshold between two windows must be between 0.0 and 1.0." << endl;
                return 1;
            }

        } else if (args[i] == "--blast_score_TH") {

            blast_score_TH = atoi(args[i + 1].c_str());

            if (blast_score_TH <= 0 || blast_score_TH > 100) {
                cerr << "BLAST hit identity and coverage score must be between 0 and 100." << endl;
                return 1;
            }

        } else if (args[i] == "--ignore_short_seq") {

            ignore_short_seq = atoi(args[i + 1].c_str());

            if (ignore_short_seq < 0) {
                cerr << "Minimum length of the contig to be removed from the assembly must be non-negative" << endl;
                return 1;
            }

        } else if (args[i] == "--insert_Ns") {

            Ns = atoi(args[i + 1].c_str());

            if (Ns < 0) {
                cerr << "Number of Ns to insert between non-merged contig and its CAP3 assembly must be non-negative" << endl;
                return 1;
            }

        } else if (args[i] == "--help") {
            help();
            exit(0);

        }

    }

    if (fastq_first == " " || fastq_sec == " " || fastq_first == fastq_sec) {

       cerr << endl << "Missing valid RNA-seq data..." << endl;
        help();
        exit(0);

    }
    
   
    if (out_directory.empty()){ // if no output directory specified by user
     size_t found = fastq_first.find_last_of("/\\");

        if (found!= string::npos) { //
               out_directory = fastq_first.substr(0, found); // set path of fastq read files for ROAST output
        } else { // running ROAST from fastq file's directory already
                out_directory = "." ;//fastq_first.substr(0, found);
        }
    
    }
   // cout << "Path is: " <<out_directory << endl;
    /*create log folder to store log files*/
    string path = out_directory + "/Log";
    char path_log[path.size() + 1];
    strcpy(path_log, path.c_str());
    filesystem::path dir(path_log);
    if (filesystem::create_directory(dir)) {
        //   std::cout << "Successfully log folder created" << "\n";	
    }
    
    path_inter = out_directory + "/Intermediate_improved_assemblies";
    char path_int[path_inter.size() + 1];
    strcpy(path_int, path_inter.c_str());
    filesystem::path dir_int(path_int);
    if (filesystem::create_directory(dir_int)) {
        //std::cout << "Successfully log folder created" << "\n";
    }


    string path_final = out_directory + "/Final_assembly";
    char path_f[path_final.size() + 1];
    strcpy(path_f, path_final.c_str());
    filesystem::path dir_final(path_f);

    if (filesystem::create_directory(dir_final)) {
        //std::cout << "Successfully log folder created" << "\n";
    }
    
    if (generate_assembly == true) {
        //develop assembly
        find = fastq_first.find_first_of(".");
        size_t find_path = fastq_first.find_last_of("/\\");
        string new_path;
        if (find_path != string::npos) {
                new_path = fastq_first.substr(0, find_path) + "/trinity_assembly";
        } else { // running ROAST from fastq file's folder
                new_path = "./trinity_assembly";
        }

        string init_fasta = fastq_first.substr(0, find_path) + "/trinity_assembly.Trinity.fasta";
        string init_ST_fasta;// = new_path + "Trinity.SuperTrans.fasta";

        //./TRINITY_HOME/Trinity --seqType fq --left SRR12002106_1.fastq --right SRR12002106_2.fastq --CPU 6 --max_memory 20G -output assembly_trinity --no_bowtie
        string trinityAssembly_command = "$TRINITY_HOME/Trinity --seqType fq --left " + fastq_first + " --right " + fastq_sec + " --normalize_reads --min_kmer_cov 5 --CPU " + threads + " --max_memory " + max_memory_TRINITY + "G --output " + new_path + " --no_bowtie  --include_supertranscripts --full_cleanup  > /dev/null 2>&1";
        // cout << trinityAssembly_command << endl;
        std::system(trinityAssembly_command.c_str());

        if (boost::filesystem::exists(init_fasta)) // does Trinity assembly file exist? generate STs
        {
            //TRINITY_HOME/Analysis/OuterTranscripts/Trinity_gene_splice_modeler.py --trinity_fasta AL.fasta --out_prefix AL.ST.fasta
            find = init_fasta.find_last_of(".");
            init_ST_fasta  = init_fasta.substr(0, find) + ".ST";
            string trinitySuperTranscript_command = "python " + exe_path + "external_tools/SuperTranscripts/Trinity_gene_splice_modeler.py --trinity_fasta " + init_fasta + " --out_prefix " + init_ST_fasta + " > /dev/null 2>&1";
           // cout << trinitySuperTranscript_command << endl;
            std::system(trinitySuperTranscript_command.c_str());

            init_ST_fasta = init_ST_fasta + ".fasta";

            if (boost::filesystem::exists(init_ST_fasta)) // does supertranscript file exist?
            {
                //cout << "continue" << endl; 
            } else {
                cerr << endl << "Couldn't generate SuperTranscripts from generated Trinity assembly. Please inspect the error or provide generated SuperTranscripts for improvement...." << endl;
                exit(0);
            }

            ST_fasta = init_ST_fasta;

            fasta = ST_fasta;
        } else {
            string trinity_version = "trinity_version";
            string trinity_version_check = "$TRINITY_HOME/Trinity --version > " + trinity_version;
            std::system(trinity_version_check.c_str());

            ifstream trinityversion;
            trinityversion.open(trinity_version.c_str()); // opens the file

            string trinity; // indata is like cin
            trinityversion >> trinity;
	    //cout << "trinity version is " << trinity << endl;
	    size_t search_trinity = trinity.find("Trinity");

            if (search_trinity != string::npos) { // Trinity file  exist
                cerr << endl << "Couldn't generate assembly properly. Please inspect Trinity dependncies error or provide generated assembly for improvement." << endl;
                exit(0);
            } else { // Trinity file doesn't found
                cerr << endl << "Couldn't generate assembly. Please set and export environmental variable TRINITY_HOME to point Trinity installation folder to generate assembly or provide generated assembly for improvement." << endl;
                exit(0);
            }
               utils.remove_file(trinity_version);
        }

    }
    if (cdhitest == 1) { // run only in the start


        //cd-hit-est -i R1.fa -o R1.95.fa  -c 0.95 -n 10 -d 0  -T 8    #To remove redundancy
        find = ST_fasta.find_last_of(".");
        ref_ST_CD_fasta = path_inter + "/initial_assembly.CD.fasta";
        CDhitEST_command = exe_path + "external_tools/cd-hit-est -i " + ST_fasta + " -o " + ref_ST_CD_fasta + "  -c 0.95  -T " + threads + " > /dev/null 2>&1";
        std::system(CDhitEST_command.c_str());

        if (boost::filesystem::exists(ref_ST_CD_fasta)) // does filePath actually exist?
        {
            //file found
        } else {

            cerr << endl << "Couldn't run cd-hit-est successfully, please check input reference file and path folder of external tools" << endl;
            cerr << "   --cdHitEST \t\t\t For any library version issue ignore removing redundancy using CD-HIT-EST by setting it to 0, default 1" << endl << endl;
            exit(0);
        }

        ifstream assem;
        assem.open(ref_ST_CD_fasta.c_str());
        if (assem.peek() == std::ifstream::traits_type::eof()) {

            cerr << endl << "Couldn't run cd-hit-est successfully, please check input reference file and path folder of external tools" << endl;
            cerr << "   --cdHitEST \t\t\t For any library version issue ignore removing redundancy using CD-HIT-EST by setting it to 0, default 1" << endl << endl;

            exit(0);
        }
        assem.close();

        fasta = ref_ST_CD_fasta;
    }


 

    /*Move fasta file name _0 to intermediate folder for further processing*/

    // command to convert multiline fasta into single line. incorporate that too in code
    //awk '{if(NR==1) {print $0} else {if($0 ~ /^>/) {print "\n"$0} else {printf $0}}}' interleaved.fasta > singleline.fasta

    string inter_assembly = path_inter + "/original_assembly.fasta";
    string command_multiTosingleLineFasta = "awk '{if(NR==1) {print $0} else {if($0 ~ /^>/) {print \"\\n\" $0} else {printf $0}}}' "  + fasta  + " > " + inter_assembly;   
    std::system(command_multiTosingleLineFasta.c_str());
    
    string init_assembly = path_inter + "/initial_assembly_0.fasta";
    string initial_IDs_map = path_inter + "/initial_IDs_map.txt";

    utils.update_contigIDs(inter_assembly, init_assembly, initial_IDs_map);
    //utils.remove_file(inter_assembly.c_str());

    find = init_assembly.find_last_of(".");

    cout << "No. of contigs in initial assembly: " << endl;
    string contig_count_command = "grep '>' " + init_assembly + " | wc -l";
    std::system(contig_count_command.c_str());

    stringstream itr;
    itr << inner_iteration;

    string fastq_first_out = path_inter + "/r1_filt.fastq" + itr.str(); //to save filtered fastq files in intermediate folder
    string fastq_sec_out = path_inter + "/r2_filt.fastq" + itr.str();

    string merged_fasta, bam_file, bam_new, cov_file, extendfasta_temp, cov_hisat, cov_minimap, merged_bam_file, inner_iteration_log, bam_sorted_filt_file;
    int merged_count, chimera_fixed_count;

    bam_new = init_assembly.substr(0, find) + ".new" + itr.str() + ".bam";
    cov_file = init_assembly.substr(0, find) + ".cov" + itr.str();



    string time_log_file = path + "/time.log";
    ofstream log_time;
    log_time.open(time_log_file.c_str());

    string byRI_log = path + "/RI.log";

    string init_fasta = init_assembly;
    string init_fastq_first = fastq_first;
    string init_fastq_sec = fastq_sec;

    bool extend = true;
    bool RI = true;
    bool improve = true;

    bySoftclip bySoftclip;

    while (outer_iteration <= outer_itr_TH && improve) { //until get improved contigs less than a TH or reached to iteration TH, keep improving

        alignment alignment;
        std::ostringstream ct;
        ct << outer_iteration;

        cout << "Outer iteration:" << outer_iteration << " started." << endl;

        while (inner_iteration <= inner_itr_TH && extend) { //&& extend


            AllFasta_data = alignment.extract_fasta_data(init_assembly); // upload fasta data in map AllFasta_data;
            contig_count = AllFasta_data.size() - 1;

            std::ostringstream itr;
            itr << inner_iteration;

            log_time << "iteration:" << inner_iteration << " started." << endl;
            cout << "Inner iteration:" << inner_iteration << " started." << endl;

            // alignment   
            time_t begin, end;
            time(&begin);

            cout << outer_iteration << "-" << inner_iteration << ": Alignment started for fragmented and partial contigs process using soft-clipping" << endl << endl;

            bam_sorted_filt_file = alignment.align_reads(init_assembly, fastq_first, fastq_sec);

            cout << "Alignment ended" << endl << endl;

            if (outer_iteration == 1 && inner_iteration == 1) {// to get average insert size only for once
                avg_IS = utils.average_insertsize(bam_sorted_filt_file);
                cout << "read length of this data is: " << read_length << endl; // " and average insert size is: " << avg_IS << endl;

                // update variables which depend on read length
                left_edge_boundary = read_length;
                right_edge_boundary = left_edge_boundary * 2;
                consecutive_missAssembled_pos_dist = read_length;
                st_end_th = read_length;
                win_size = read_length * 2;                
            }

            time(&end);
            double elapsed_secs = double(end - begin);
            log_time << "Outer iteration:" << outer_iteration << " Inner iteration:" << inner_iteration << " alignment time:" << elapsed_secs << endl;

            //extend assembly
            time_t begin2, end2;
            time(&begin2);

            cout << outer_iteration << "-" << inner_iteration << ": Assembly extension using soft clipped reads started" << endl;

            merged_fasta = path_inter + "/Improved_assembly_SCs_" + ct.str() + "-" + itr.str() + ".fasta";
            inner_iteration_log = path + "/iteration_" + ct.str() + "_" + itr.str();

            extend = bySoftclip.extend_bySoftclip(init_assembly, merged_fasta, bam_sorted_filt_file, inner_iteration_log);

            cout << "Assembly extension using soft clipped reads ended" << endl << endl;

            time(&end2);
            double elapsed_secs2 = double(end2 - begin2);
            log_time << "Outer iteration:" << outer_iteration << " Inner iteration:" << inner_iteration << " merging &  extension time:" << elapsed_secs2 << endl;

            if (cdhitest == 2) {
                //cd-hit-est -i R1.fa -o R1.95.fa  -c 0.95 -n 10 -d 0  -T 8    #To remove redundancy
                find = merged_fasta.find_last_of(".");
                temp_assembly = merged_fasta.substr(0, find) + ".CD.fasta";
                CDhitEST_command = exe_path + "external_tools/cd-hit-est -i " + merged_fasta + " -o " + temp_assembly + "  -c 0.95  -T " + threads + " > /dev/null 2>&1";
                std::system(CDhitEST_command.c_str());

                // init_assembly = temp_assembly; //merged_fasta;

                string duplication_remove = "awk '/^>/{f=!d[$1];d[$1]=1}f' " + temp_assembly + " > " + init_assembly;
                std::system(duplication_remove.c_str());
                //utils.remove_file(merged_fasta.c_str());
            } else {
                string duplication_remove = "awk '/^>/{f=!d[$1];d[$1]=1}f' " + merged_fasta + " > " + init_assembly; //temp_assembly + " > " + init_assembly;
                std::system(duplication_remove.c_str());
            }
                utils.remove_file(temp_assembly.c_str());
            // merged_fasta = init_fasta.substr(0, init_fasta.size()-5) + "ext_merg.fasta"+ itr.str();

            //    if(extend){
            //filter sam -> fastq
            cout << outer_iteration << "-" << inner_iteration << ": Fastq files filtering process started." << endl;
            time_t begin3, end3;
            time(&begin3);
            filterSamToFastq fil_sam_to_fastq;
            fil_sam_to_fastq.filterSam(bam_sorted_filt_file, bam_new, fastq_first, fastq_first_out, fastq_sec, fastq_sec_out);
            cout << "Fastq file filtering step ended" << endl << endl;


            time(&end3);
            double elapsed_secs3 = double(end3 - begin3);
            log_time << "Outer iteration:" << outer_iteration << " Inner iteration:" << inner_iteration << " fastq filter time:" << elapsed_secs3 << endl << endl;

            //// fil_sam_to_fastq.remove_prev_fastq(fastq_first, fastq_sec);
            fastq_first = fastq_first_out;
            fastq_sec = fastq_sec_out;

            fastq_first_out = path_inter + "/r1_filt.fastq" + itr.str(); //to save filtered fastq files in intermediate folder
            fastq_sec_out = path_inter + "/r2_filt.fastq" + itr.str();

            utils.remove_file(cov_file.c_str());

            // log_time << "Go back to alignment step" << endl << endl ;
            std::ifstream fastq1, fastq2;
            string entry;
            fastq1.open(fastq_first.c_str());
            fastq2.open(fastq_sec.c_str());

            if (getline(fastq1, entry).eof() || getline(fastq2, entry).eof()) // add check when filtered bam file is too small to produce fastq file out of it 
            {
                fastq_first = init_fastq_first;
                fastq_sec = init_fastq_sec;
            }
            fastq1.close(); fastq2.close();
            inner_iteration++;
            //}
            //      else
            //       goto softclip_process_finish;
            //filter
            //AllFasta_data.clear();
        }
        //softclip_process_finish:
        //     ;

        utils.remove_file(cov_file.c_str());
        //filterSamToFastq fil_sam_to_fastq;
        //fil_sam_to_fastq.remove_prev_fastq(fastq_first, fastq_sec);

        /*fix fragmented contigs again which are extended by softclips */

        string mergedSCsfasta = path_inter + "/Improved_assembly_SCs_merged_" + ct.str() + ".fasta";
        string sc_merged_log = path + "/SCs_merged_iteration_" + ct.str() + ".log";

        cout << "***** Assembly merging using soft-clipped-extended ends: Process started *****" << endl << endl;

        int merged_count_SCs = bySoftclip.fix_fragmentedBySCs(init_assembly, mergedSCsfasta, sc_merged_log);
        //.FindFragmentsBySCs("1", assembly, path_inter, exe_path, SC_cons_len_TH, min_overlap_TH, BLAST_score_TH, max_allowed_gaps);

        //exit(0);
        string delimiter = "\t";

        cout << "***** Assembly extension and merging using soft clipped reads: Process completed *****" << endl << endl;

        if (RI) {
            utils.remove_file(bam_new.c_str());
            //}
            /*correct assembly using unmapped and distantly mapped reads*/
            time_t begin4, end4;
            time(&begin4);
            cout << outer_iteration << ": Re-alignment on updated assembly for fragmented and partial contigs process using Unmapped and discordantly mapped reads" << endl << endl;
            merged_fasta = mergedSCsfasta;

            AllFasta_data = alignment.extract_fasta_data(merged_fasta); // upload fasta data in map AllFasta_data;
            contig_count = AllFasta_data.size() - 1;
            // cout << "NO. of contigs" << contig_count <<endl;

            string new_bam_sorted_filt_file;
            //  if (outer_iteration == 1)
            // new_bam_sorted_filt_file = "/media/madiha/Data3/DTR/Rice/rice_benchmark/intermediate_Improved_assemblies/initial_assembly_0.sorted.filtered.bam";
            // else
            new_bam_sorted_filt_file = alignment.align_reads(merged_fasta, fastq_first, fastq_sec); // cov_file); 
            //call byRI with new fasta, bam and log file

            time(&end4);
            double elapsed_secs4 = double(end4 - begin4);
            log_time << "Outer iteration:" << outer_iteration << " Unmapped and discordantly mapped reads process alignment time:" << elapsed_secs4 << endl;

            cout << outer_iteration << ": Assembly extension and merging using Unmapped and discordantly mapped reads: process started" << endl;
            byRI ri;

            //  time_t begin5, end5;
            // time(&begin5);
            extendfasta_temp = path_inter + "/Improved_assembly_RI_temp" + ".fasta";

            merged_count = ri.read_island(new_bam_sorted_filt_file, merged_fasta, extendfasta_temp, byRI_log, log_time);
            // utils.remove_file(new_bam_sorted_filt_file);
            //utils.remove_file(cov_file);

            cout << "*****  Assembly extension and merging using Unmapped and discordantly mapped reads: process completed  *****" << endl << endl;
            // AllFasta_data.clear();
            // time(&end5);
            // double elapsed_secs5 = double(end5 - begin5);
            // log_time << "outer_iteration:" << outer_iteration << " Read islands process time:" << elapsed_secs5 << endl << endl;

        } else
            extendfasta_temp = init_assembly;
        time_t begin6, end6;
        time(&begin6);

        string extendfasta = path_inter + "/Improved_assembly_RI_" + ct.str() + ".fasta";


        //cd-hit-est -i R1.fa -o R1.95.fa  -c 0.95 -n 10 -d 0  -T 8    #To remove redundancy
        find = extendfasta_temp.find_last_of(".");
        
        if (cdhitest == 2) {
            temp_assembly = extendfasta_temp.substr(0, find) + ".CD.fasta";
            CDhitEST_command = exe_path + "external_tools/cd-hit-est -i " + extendfasta_temp + " -o " + temp_assembly + "  -c 0.95  -T " + threads + " > /dev/null 2>&1";
            std::system(CDhitEST_command.c_str());

            //remove duplications from output
            //awk '/^>/{f=!d[$1];d[$1]=1}f' Improved_assembly_RI_4.fasta > out.fa
            string duplication_remove = "awk '/^>/{f=!d[$1];d[$1]=1}f' " + temp_assembly + " > " + extendfasta;
            std::system(duplication_remove.c_str());
        } else {
            string duplication_remove = "awk '/^>/{f=!d[$1];d[$1]=1}f' " + extendfasta_temp + " > " + extendfasta; //temp_assembly + " > " + extendfasta;
            std::system(duplication_remove.c_str());
        }
        utils.remove_file(extendfasta_temp.c_str());
        utils.remove_file(temp_assembly.c_str());

        /*correct local mis-assemblies and false chimeras*/
        cout << outer_iteration << ": Re-alignment on updated assembly for Local mis-assemblies/false chimeras identification process" << endl << endl;

        AllFasta_data = alignment.extract_fasta_data(extendfasta); // upload fasta data in map AllFasta_data;
        contig_count = AllFasta_data.size() - 1;

        find = extendfasta.find_last_of(".");
        cov_minimap = extendfasta.substr(0, find) + ".minimap2.cov" + itr.str();

        cov_hisat = extendfasta.substr(0, find) + ".hisat.cov" + itr.str();

        merged_bam_file = alignment.align_reads(extendfasta, init_fastq_first, init_fastq_sec, cov_minimap, cov_hisat); // minimap2 alignment return only because from hisat2 we need coverage only


        time(&end6);
        double elapsed_secs6 = double(end6 - begin6);
        log_time << "Outer iteration:" << outer_iteration << " Local mis-assemblies/false chimeras process alignment time:" << elapsed_secs6 << endl;


        time_t begin7, end7;
        time(&begin7);

        cout << outer_iteration << ": Identification of Local mis-assemblies/false chimeras: Process started" << endl;
        string mis_assembly_chimera_log = path + "/mis_assembly_chimera" + ct.str() + ".log";
        string processed_fastafile = path_inter + "/Improved_assembly_MC_" + ct.str() + ".fasta";
        string cufflink = path_inter + "/transcripts.gtf";


        mis_assembly_chimera chimera;
        chimera_fixed_count = chimera.extract_chimera_mis_assembly(merged_bam_file, cov_hisat, cov_minimap, extendfasta, processed_fastafile, cufflink, mis_assembly_chimera_log);

        cout << "*****  Identification of Local mis-assemblies/false chimeras: Process completed  *****" << endl << endl;


        time(&end7);
        double elapsed_secs7 = double(end7 - begin7);
        log_time << "Outer iteration:" << outer_iteration << " Identification of Local mis-assemblies/false chimeras process time:" << elapsed_secs7 << endl << endl;


        utils.remove_file(merged_bam_file);
        utils.remove_file(cov_hisat);
        utils.remove_file(cov_minimap);
        utils.remove_file(cufflink);

        //cd-hit-est -i R1.fa -o R1.95.fa  -c 0.95 -n 10 -d 0  -T 8    #To remove redundancy
        find = processed_fastafile.find_last_of(".");

        if (cdhitest == 2) {
            temp_assembly = processed_fastafile.substr(0, find) + ".CD.fasta";
            CDhitEST_command = exe_path + "external_tools/cd-hit-est -i " + processed_fastafile + " -o " + temp_assembly + "  -c 0.95  -T " + threads + " > /dev/null 2>&1";
            std::system(CDhitEST_command.c_str());

            //remove duplications from output
            //awk '/^>/{f=!d[$1];d[$1]=1}f' Improved_assembly_RI_4.fasta > out.fa
            string duplication_remove = "awk '/^>/{f=!d[$1];d[$1]=1}f' " + temp_assembly + " > " + init_assembly;
            std::system(duplication_remove.c_str());

        string remove_CDhitfiles = "rm " + path_inter + "/*.clstr";
        std::system(remove_CDhitfiles.c_str());
        } else {
            string duplication_remove = "awk '/^>/{f=!d[$1];d[$1]=1}f' " + processed_fastafile + " > " + init_assembly; //temp_assembly + " > " + init_assembly;
            std::system(duplication_remove.c_str());
        }
        utils.remove_file(temp_assembly.c_str());


        // init_assembly = processed_fastafile;
        fastq_first = init_fastq_first;
        fastq_sec = init_fastq_sec;
        inner_iteration = 1;
        extend = true;
        RI = true;
        outer_iteration++;

        if (merged_count >= improvment_TH || merged_count_SCs >= improvment_TH || chimera_fixed_count >= improvment_TH) //untill get improved contigs less than a TH, keep improving
            improve = true;
        else
            improve = false;

        //AllFasta_data.clear();
        //delete everything except recently updated assembly

        string remove_allBAM = "rm " + path_inter + "/*.bam";
        std::system(remove_allBAM.c_str());

        string remove_allBAMindex = "rm " + path_inter + "/*.bai";
        std::system(remove_allBAMindex.c_str());

        string remove_allfilteredfastq = "rm " + path_inter + "/*.fastq*";
        std::system(remove_allfilteredfastq.c_str());

        /* stop file found in the folder? terminate this iteration and go out of outer loop*/

        string stop_ROAST = path_inter + "/stop.txt";

        if (boost::filesystem::exists(stop_ROAST)) // does stop file exist?
        {
            break; // go out of while loop
        }

    }
    log_time.close();

    cout << endl << endl;
    cout << "*****************************************************************" << endl << "                         Process completed                         " << endl << "*****************************************************************" << endl << endl;

    //format headers of final improved fasta files
    cout << "No. of contigs in initial assembly: " << endl;
    contig_count_command = "grep '>' " + fasta + " | wc -l";
    std::system(contig_count_command.c_str());

    cout << "Total number of merges occurred: " << merged_fragmented_contigs << endl;
    cout << "Total number of splits occurred: " << fixed_chimeric_contigs << endl;

    //format headers of final improved fasta files
    string final_fasta = path_final + "/final_assembly.fasta";
    utils.format_finalAssembly(init_assembly, final_fasta, initial_IDs_map);

    cout << "No. of contigs in Final assembly: " << endl;
    contig_count_command = "grep '>' " + final_fasta + " | wc -l";
    std::system(contig_count_command.c_str());

    cout << "Total number of extensions occurred: " << extended_incomplete_contigs << endl;
    cout << "Total number of updated sequences: " << updatad_local_misAssemblies << endl << endl;

    cout << "*****************************************************************" << endl;
    cout << "*****************************************************************" << endl;


    if (complete_cleanup) {
        boost::filesystem::remove_all(path);
        boost::filesystem::remove_all(path_inter);
    }
    return 0;
}