%TSS and PAS are retrieved from known genes
%zheng xia@tmhs.org, Dec-26th-2010

clear;clc;
global epsilon;
epsilon = 1e-10;

%input file name
input_filename = 'test.instances';
output_filename = 'test_result';


fid = fopen(input_filename);
fid_write = fopen(output_filename, 'wt');
NN = fscanf(fid,'%f', 1);       
read_len = fscanf(fid,'%f', 1);  

while(1)

    name = fscanf(fid,'%s',1);      
    
    if (strcmp(name,''))         
        break;
    end
    m = fscanf(fid,'%d',1);      
    
    %read the annotated isoforms
    n = fscanf(fid,'%d',1);      %isoform number of annotated isoforms
    annotated_isoforms = cell(n,5);
    for ni = 1:n
        for nj=1:5                    
            annotated_isoforms{ni,nj} = fscanf(fid,'%s',1); %isoform name
        end
    end
    
   
    
    
    
    tophat_exons = fscanf(fid,'%d',2*m); %positions from tophat
    tophat_exons = reshape(tophat_exons,m,2); 
    
    %TSS_PASs
    temp_label = fscanf(fid,'%s',1);
    num_TSSPASs = fscanf(fid,'%d',1);
    TSS_PASs = zeros(num_TSSPASs,2);
    for ni = 1:num_TSSPASs
        TSS_PASs(ni,1) = fscanf(fid,'%d',1)+1;
        TSS_PASs(ni,2) = fscanf(fid,'%d',1)+1;
    end
    
    
    
    L = fscanf(fid,'%f',m)';           %exon length
    N = fscanf(fid,'%f',m)';           %number of reads on each exon 
    if m>1
        L1 = fscanf(fid,'%f',m*(m-1)/2)'; %alternative splice length
        N1 = fscanf(fid,'%f',m*(m-1)/2)'; %number of reads on the junction region
    else
        L1 = [];
        N1 = [];
    end


    
    if m > 1 && m < 80   
        
        name
     
        p = 0.5;  %sparsity 
        [results, A_candidate,list_junction] = BLasso_whole_path_real_shortestpath_TSSPAS(TSS_PASs,tophat_exons,L,L1,N,N1,NN,...
                                                                               read_len,name,p);
                                                                         
                                                                           
        betas2 = results.L1_nbetas; 

        ites = size(betas2,2);
        Inds_selected = find(betas2(:,1)>0); 
        max_isoform_num = 8;  %The maximum number for expressed isoforms

        for i=2:ites
            temp = find(betas2(:,i)>0);
            Inds_selected = setcheck(temp,Inds_selected,max_isoform_num);
        end

        m_inds = length(Inds_selected);
        objs_val = zeros(1,m_inds);
        model_len = zeros(1,m_inds); 
        for i = 1:m_inds
            [X_ests, objs_val(i)] = isoform_estimation_mine(A_candidate(:,Inds_selected{i}),L,L1,N,N1,NN,read_len,name);
            model_len(i) = length(Inds_selected{i});
        end


        model_len_unique = unique(model_len); 

        objs_val2 = []; 
        length_model_candid = length(model_len_unique);
        model_candidate = [];
        for i = 1:length_model_candid,
            [temp temp_ind] = max(objs_val(model_len==model_len_unique(i)));
            temp_inds = find(model_len==model_len_unique(i));
            model_candidate = [model_candidate temp_inds(temp_ind)];
            objs_val2 = [objs_val2 temp];

        end

        %%%%%%%%%%%%%%%%%%%
        %ratio test
        model_length = length(model_candidate);
        model_sel = 1;
        for i = 2:model_length
            H_test = lratiotest(objs_val2(i),objs_val2(i-1),1,1e-20);
            if H_test ~= 1
                model_sel = i-1;
                i = i - 1;
                break;
            end
        end

        if i == model_length   
            model_sel = model_length;
        end

        if isempty(model_candidate)
           X_selected = 0; 
           ind_sel2   = [];
        else
            model_sel2 = model_candidate(model_sel);
            ind_sel2 = Inds_selected{model_sel2};
            [X_selected, b] = isoform_estimation_mine(A_candidate(:,ind_sel2),L,L1,N,N1,NN,read_len,name);
        end
    
        
        %%%%%%%%%%%%%%%%%%%%%%%
        if size(ind_sel2,1)>0
            isoforms_result = A_candidate(:,ind_sel2);
            chromosome = annotated_isoforms{2};
            strand = annotated_isoforms{3};
            

            for i = 1:length(ind_sel2)
                isoform_exons = isoforms_result(:,i);
                isoform_exon_ind = find(isoform_exons==1);
                transcript_start = num2str(tophat_exons(isoform_exon_ind(1),1)+1);
                transcript_end   = num2str(tophat_exons(isoform_exon_ind(end),2));
                %write the transript
                fprintf(fid_write, '%s\tNSMAP\ttranscript\t%s\t%s\t1000\t%s\t.\tgene_id \"NSMAP.%s\"; transcript_id \"NSMAP.%s.%d\"; RPKM \"%g\"\n',...
                         chromosome,transcript_start,transcript_end,strand,name,name,i,X_selected(i));
                
                %write the exon information
                exons_selected_in_cur_isoform = tophat_exons(isoform_exon_ind,:);
                
                exons_processed = exons_selected_in_cur_isoform(1,:);
                for j = 2: length(isoform_exon_ind)
                    if exons_selected_in_cur_isoform(j,1) - exons_processed(end,2) == 0
                        exons_processed(end,2) = exons_selected_in_cur_isoform(j,2);
                    else
                        exons_processed = [exons_processed;exons_selected_in_cur_isoform(j,:)];
                    end
                end
                
                nums_exons_cur_isoform = size(exons_processed,1);
                for j = 1:nums_exons_cur_isoform
                    fprintf(fid_write, '%s\tNSMAP\texon\t%s\t%s\t1000\t%s\t.\tgene_id \"NSMAP.%s\"; transcript_id \"NSMAP.%s.%d\"\n',...
                         chromosome,num2str(exons_processed(j,1)+1),num2str(exons_processed(j,2)),strand,name,name,i);
                end
                     
                
            end
        end
        
     elseif m==1

         cur_isoform_RPKM = N*1e9/NN/L;
         chromosome = annotated_isoforms{2};
         strand = annotated_isoforms{3};
         transcript_start = num2str(tophat_exons(1,1)+1);
         transcript_end   = num2str(tophat_exons(1,2));
         %write the transript
         fprintf(fid_write, '%s\tNSMAP\ttranscript\t%s\t%s\t1000\t%s\t.\tgene_id \"NSMAP.%s\"; transcript_id \"NSMAP.%s.%d\"; RPKM \"%g\"\n',...
                         chromosome,transcript_start,transcript_end,strand,name,name,1,cur_isoform_RPKM);
         fprintf(fid_write, '%s\tNSMAP\texon\t%s\t%s\t1000\t%s\t.\tgene_id \"NSMAP.%s\"; transcript_id \"NSMAP.%s.%d\"\n',...
                         chromosome,transcript_start,transcript_end,strand,name,name,1);
        
        
     end

end
fclose(fid);
fclose(fid_write);

