%% Explicació

%Aquest script agafa cadascuna de les proteïnes trobades en Treball Bioinfo
%_3 > Proteins_musmusculus i hi crea una figura de matlab per a cadascun
%dels contigs en els que BLAST ha trobat algun hit.

%Aquesta figura ensenya el contig, els hits del blast (en cadena positiva i
%negativa), els genms predits per exonerate/genewise i la part dels exons a
%la qual fa referencia cada part predita i els hits del blast


%% Script que ho fa tot de cop

clear all

%Defineix el llistat de Selenoproteïnes de Mus Musculus (guardades en All_folders)

%%%%%%%%% Cada PC és diferent
project_dir = [ getenv('C:') filesep 'Users' filesep 'mikis_000' filesep 'Desktop' filesep ] ; %Miki PC Exemple

%%%%%%%%%%

%Llistat de carpetes on hi ha les proteïnes:
All_files_DS = struct2ds(dir(project_dir))
All_folders = All_files_DS(find(regexpcmp(All_files_DS.name,'SPP') == 1),:).name;

% Entra a cada carpeta
for Folder = All_folders' %cada cop entra en una carpeta
% for Folder = {'SELENOV'} %Carpeta Exemple 
    Folder
    
    if isempty(find(regexpcmp(Folder{1},'_'))) == 1
    %%%%%%%%%
    direct = [strcat(project_dir ,Folder{1})]; %Miki PC
    %%%%%%%%%
    
    cd(direct); %posa't al directori de la carpeta
    Content = struct2ds(dir(direct));
    Content = Content.name; %nom de tots els fitxers en la carpeta
    
    if length(Content) > 4
    Blast = importdata('filtrat_blast.fa');
    Blast_Starts = Blast.data(:,1);
    Blast_Ends = Blast.data(:,2);
    Blast_Evals = Blast.data(:,3);
    All_Contigs = Blast.textdata; %contigs trobats en el blast
    Unique_Contigs = unique(All_Contigs);
    
    %una figura amb tots els contigs:
    
    SumToY = 0; %volem posar tots els contigs en la mateixa figura, i això vold dir que cada cop anirem cap amunt en l'eixY 
    fig = figure; hold on; grid on;    
    
    for Contig = Unique_Contigs' 
    
%     for Contig = Unique_Contigs(7) 
%         Contig
        idx_cont = find(regexpcmp(Content,['subseq_',Contig{1}]));
        if isempty(idx_cont) == 0 %si els hits d'aquest contig han propsperat(E<0.01)
%         Cont_Info = fastainfo(Content{idx_cont(1)});
%         C_Len = double(Cont_Info.Length);
%         C_Len = 500000000; %if fasta length fails
        idx_len = find(regexpcmp(Content,['length_',Contig{1}]));
        
        if isempty(idx_len) == 0
        C_Len_File = importdata(Content{idx_len(1)});
        C_Len_File = C_Len_File{1};
        Error = 1;
        try
        SpaceID = regexp(C_Len_File,'C');
        C_Len = str2num(C_Len_File(1:SpaceID-2));
        linspace(1,C_Len,4);
        Error = 0;
        end
        if Error == 1
        SpaceID = regexp(C_Len_File,'L');    
        end
        C_Len = str2num(C_Len_File(1:SpaceID-2));

        plot(linspace(1,C_Len,4),linspace(SumToY,SumToY,4),'color','k'); %posem la linia del contig a 0
        text(0,0.2+SumToY,Contig)
        % posa els hits del blast:
        idx_blast = find(strcmp(All_Contigs,Contig)); %defineix l'índex de  de cadascun dels contigs per fer els hits
        for B_h = idx_blast' % ves per tots els hits del blast
            First = Blast_Starts(B_h);
            Second = Blast_Ends(B_h);
            if First < Second %si la orientació és FWD
                Pos = 0.1 + SumToY;
                Start = First;
                End = Second; 
            end 
            if First > Second %si la orientació és FWD
                Pos = -0.2 + SumToY;
                Start = Second;
                End = First;   
            end            
            rectangle('Position',[Start  Pos End-Start 0.1],'FaceColor','m','EdgeColor','m')%fes un rectangle per cada hit [x,y,w,h]
        end
        %posa els subseq predit a cada contig
        idx_subseqs = find(regexpcmp(Content,['subseq_',Contig{1}])); 
        ss_starts = []; %vector que compila tots els starts dels subseqs
        ss_ends = [];
        S = 1;
        for ss = Content(idx_subseqs)' %ves per cada subseq
            ss_file = importdata(ss{1});
            File_Name = ss{1};
            ss_ID = str2num(File_Name(length(File_Name)-3));
            Name = ss_file{1,1};
            pos_par_open = regexp(Name,'('); %pos del(
            pos_par_close = regexp(Name,')');
            pos_coma = regexp(Name,',');
            ss_Start = str2num(Name(pos_par_open+1:pos_coma-1));
            ss_End = ss_Start + str2num(Name(pos_coma+1:pos_par_close-1));
            rectangle('Position',[ss_Start  -0.1+SumToY ss_End-ss_Start 0.2],'FaceColor','k','EdgeColor','k')
            text(ss_Start-100,SumToY,strcat('ssID =',{' '},num2str(ss_ID)))
            S = S + 1;   
            ss_starts = [ss_starts ; ss_Start];
            ss_ends = [ss_ends ; ss_End];
            %posa els gens (separats per exons) predits per exonerate (mira
            %que siguin del mateix ID):
            idx_exonerates = find(regexpcmp(Content,['exonerate_final_',Contig{1}]));
            for gene_file = Content(idx_exonerates)'
                gene_filename = gene_file{1};
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CANVIAR AMB NOU SET (el 6 serà 5)%%              
                ss_gene_ID = str2num(gene_filename(length(gene_filename)-5));
                gene_ID = str2num(gene_filename(length(gene_filename)-4));
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
                Gene_Struct = GFFAnnotation(gene_file{1});
                Starts = double(Gene_Struct.Start);
                Ends = double(Gene_Struct.Stop);
                Strand = unique(Gene_Struct.Strand);
                Attributes = Gene_Struct.Attributes;
                if strcmp(Strand,'+') == 1; Pos = 0.2 + SumToY; sumMid = 0.05; end
                if strcmp(Strand,'-') == 1; Pos = -0.3 + SumToY; sumMid = 0.05; end
                
                if ss_gene_ID == ss_ID & isempty(Starts) == 0% si es tracta d'un gen predit en aquella subseq
                    for E = 1:length(Starts) %per cada exó
                        e_Start = ss_Start + Starts(E);
                        e_length = Ends(E) - Starts(E);
%                         At = Attributes(E)
%                         text(e_Start,Pos+0.1,At)                        
                        rectangle('Position',[e_Start  Pos e_length 0.1],'FaceColor','b','EdgeColor','b')
                        if E == 1   
                        end                      
                    end
                    G_Start = min(ss_Start+Starts);
                    G_End = max(ss_Start+Ends);
                    text(G_Start,Pos+0.07,strcat('Ex-gID =',{' '},num2str(gene_ID)))
                    text(G_Start,Pos+0.03,strcat('Coord:',{' '},num2str(G_Start),'-',num2str(G_End)))
                    for In = 1:length(Starts) - 1 %per cada intró
                        i_Start = ss_Start + Ends(In)+ 1;
                        i_End = ss_Start + Starts(In+1)- 1;   
                        plot([i_Start;i_End],[Pos+sumMid;Pos+sumMid],'-','color','b','linewidth',10)            
                    end
                    
                    %afegeix t-coffee de l'exonerate:
                    %anirà entre [0.4,0.5 o -0.4,-0.5]
                    idx_tc = find(regexpcmp(Content,['tcoffee_exonerate_',Contig{1},'_',num2str(ss_ID),num2str(gene_ID)]));
                    filename = Content{idx_tc};
                    FileDS = dataset('file',filename,'Delimiter','\t');
                    TC_content = dataset2cell(FileDS);
                    if strcmp(Strand,'+') == 1
                        Pos = 0.4 + SumToY;  
                        Sum_aa = 3;
                        tc_First = G_Start;

                    end
                    if strcmp(Strand,'-') == 1 
                        Pos = -0.5 + SumToY;
                        Sum_aa = -3;
                        tc_First = G_End-3;
                    end
                    aa_Start = tc_First;
                    LenProt = 0; %llargada de la proteïna de Mus Musculus
                    N_matches = 0;
                    SelCysMatch = 'No SelCys';
                    for tc = 1:length(TC_content) %ves per totes les files del tcoffe
                        TC = TC_content{tc};
                        
                        if length(TC) > 3
                        if strcmp(TC(1:3),'SPP') == 1  %si trobes la seq de la proteïna en aquella fila
                            TC_pred = TC_content{tc+1};
                            start_prot = max(regexp(TC,'\s')) + 1;
                            start_pred = max(regexp(TC,'\s')) + 1;
                            pept_prot = TC(start_prot:length(TC));
                            pept_pred = TC_pred(start_pred:length(TC_pred));   
                            LenProt = LenProt + length(pept_prot);
                            for P = 1:length(pept_prot) %per cada aa fas un rectangle
                                aa_prot = pept_prot(P);
                                aa_pred = pept_pred(P);
                                if strcmp(aa_prot,aa_pred) == 1; color = 'g'; N_matches = N_matches + 1; end %match
                                if strcmp(aa_prot,aa_pred) == 0; color = 'r'; end %miss or in/del
                                if strcmp(aa_pred,'-') == 1; color = 'y'; end %deletion
                                if strcmp(aa_prot,'-') == 1; color = 'm'; end %insertion
                                rectangle('Position',[aa_Start  Pos 3 0.1],'FaceColor',color,'EdgeColor',color)
                                % add Se-Cys properties
                                if strcmp(aa_prot,'X') == 1
                                  r_Se_Cys = rectangle('Position',[aa_Start  Pos 3 0.1],'FaceColor',color,'EdgeColor','b');
                                  uistack(r_Se_Cys,'top')
                                  if strcmp(aa_prot,aa_pred) == 1
                                      SelCysMatch = 'true';

                                  else
                                      SelCysMatch = [num2str(aa_prot),'>',num2str(aa_pred)];
                                  end
                                if strcmp(aa_pred,'-') == 1
                                     SelCysMatch = 'true';
                                end
                                end
                                aa_Start = aa_Start + Sum_aa;
                            end      
                        end
                        end
                    end
                    %add t_coffee score
                    TC_Data = importdata(filename);
                    Info = TC_Data{1};
                    Score_start = regexp(Info,'SCORE=')+6;
                    RemainingInfo = Info(Score_start:length(Info));
                    AllComaPos = regexp(RemainingInfo,',');
                    ComaPos = AllComaPos(1);
                    Score = RemainingInfo(1:ComaPos-1);
                    % %of matches
                    pct_match_num = (N_matches/LenProt)*100;
                    pct_match = num2str(pct_match_num);
                    t_ex = text(tc_First,Pos+0.05,strcat('Ex-TC-Score:',Score,{' '},pct_match,'% id',{' '},'Al-SeCys:',SelCysMatch));
                    
                    if pct_match_num <= 30
                        set(t_ex,'color','r')
                    end     
                    if pct_match_num > 30 & pct_match_num <= 60
                        set(t_ex,'color','m')
                    end  
                    if pct_match_num > 60 & pct_match_num < 90
                        set(t_ex,'color','b')
                    end                      
                    if pct_match_num >= 90
                        set(t_ex,'color','g')
                    end
                    
                    
                   %%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%
                end                  
            end   

            %posa tot lo que té a veure amb GeneWise:
            %posa els gens
            idx_GW = find(regexpcmp(Content,['genewise_',Contig{1}]));

            for file = Content(idx_GW)
                filename = file{1};  
                ss_gene_ID = str2num(filename(length(filename)-3));                
                if ss_gene_ID == ss_ID

                    GW_output = dataset('file',filename,'Delimiter','\t');
                    if double(isempty(GW_output)) == 0
                    idx_genes = (find(strcmp(GW_output.genewiseOutput,'match')));
                    Ngenes = length(idx_genes);
                    
                    for g = 1:Ngenes %per cadascun dels gens predits per GW
                        if g+1 <= Ngenes
                            Gene_Feats = GW_output.genewiseOutput(idx_genes(g):idx_genes(g+1)-1);
                        else
                            Gene_Feats = GW_output.genewiseOutput(idx_genes(g):length(GW_output.genewiseOutput));
                        end                        
                        G_Coords = [str2num(Gene_Feats{2});str2num(Gene_Feats{3})];
                        G_Start = ss_Start + min(G_Coords);
                        G_End = ss_Start + max(G_Coords);               
                        Strand = Gene_Feats(5);
                        if strcmp(Strand,'+') == 1; Pos = 0.3 + SumToY; sumMid = 0.05; end
                        if strcmp(Strand,'-') == 1; Pos = -0.4 + SumToY; sumMid = 0.05; end
                        for F = 1:length(Gene_Feats) %ves per totes les columnes referents a aquell gen
                            if strcmp(Gene_Feats(F),'cds') == 1 %si troba un exó
                                e_Coords = [str2num(Gene_Feats{F+1});str2num(Gene_Feats{F+2})];
                                e_Start = ss_Start + min(e_Coords);
                                e_End = ss_Start + max(e_Coords);
                                e_length = e_End - e_Start;
                                rectangle('Position',[e_Start  Pos e_length 0.1],'FaceColor','y','EdgeColor','y')
                            end
                            if strcmp(Gene_Feats(F),'intron') == 1 %si troba un intró
                                i_Coords = [str2num(Gene_Feats{F+1});str2num(Gene_Feats{F+2})];
                                i_Start = ss_Start + min(i_Coords);
                                i_End = ss_Start + max(i_Coords);
                                plot([i_Start;i_End],[Pos+sumMid;Pos+sumMid],'-','color','y','linewidth',10)
                            end
                        end
                        text(G_Start,Pos+sumMid,strcat('Coord:',{' '},num2str(G_Start),'-',num2str(G_End)));

                        %tcoffee del GeneWise per cada gen trobat
                        idx_tc = find(regexpcmp(Content,['t_coffee_',Contig{1}]) & regexpcmp(Content,['_',num2str(ss_ID)]));

                        if double(isempty(idx_tc)) == 0
                        filename = Content{idx_tc};

                        FileDS = dataset('file',filename,'Delimiter','\t');
                        TC_content = dataset2cell(FileDS);
                        if strcmp(Strand,'+') == 1
                            Pos = 0.5 + SumToY;  
                            Sum_aa = 3;
                            tc_First = G_Start;
                            
                        end
                        if strcmp(Strand,'-') == 1 
                            Pos = -0.6 + SumToY;
                            Sum_aa = -3;
                            tc_First = G_End-3;
                        end

                        aa_Start = tc_First;
                        LenProt = 0; %llargada de la proteïna de Mus Musculus
                        N_matches = 0;
                        SelCysMatch = 'No SelCys';
                        for tc = 1:length(TC_content) %ves per totes les files del tcoffe
                            TC = TC_content{tc};
                            if length(TC) > 3
                            if strcmp(TC(1:3),'SPP') == 1 %si trobes la seq de la proteïna en aquella fila
                                TC_pred = TC_content{tc+1};
                                start_prot = max(regexp(TC,'\s')) + 1;
                                start_pred = max(regexp(TC,'\s')) + 1;
                                pept_prot = TC(start_prot:length(TC));
                                pept_pred = TC_pred(start_pred:length(TC_pred));          
                                LenProt = LenProt + length(pept_prot);
                                for P = 1:length(pept_prot) %per cada aa fas un rectangle
                                    aa_prot = pept_prot(P);
                                    aa_pred = pept_pred(P);
                                    if strcmp(aa_prot,aa_pred) == 1; color = 'g'; N_matches = N_matches + 1; end %match
                                    if strcmp(aa_prot,aa_pred) == 0; color = 'r'; end %miss or in/del
                                    if strcmp(aa_pred,'-') == 1; color = 'y'; end %deletion
                                    if strcmp(aa_prot,'-') == 1; color = 'm'; end %insertion
                                    rectangle('Position',[aa_Start  Pos 3 0.1],'FaceColor',color,'EdgeColor',color)
                                    % add Se-Cys properties
                                    if strcmp(aa_prot,'X') == 1
                                      r_Se_Cys = rectangle('Position',[aa_Start  Pos 3 0.1],'FaceColor',color,'EdgeColor','b');
                                      uistack(r_Se_Cys,'top')
                                      if strcmp(aa_prot,aa_pred) == 1
                                          SelCysMatch = 'true';
                                      else
                                          SelCysMatch = [num2str(aa_prot),'>',num2str(aa_pred)];
                                      end
                                    end
                                    aa_Start = aa_Start + Sum_aa;
                                end      
                            end 
                            end
                        end
                    
                        %add t_coffee score
                        TC_Data = importdata(filename);
                        Info = TC_Data{1};
                        Score_start = regexp(Info,'SCORE=')+6;
                        RemainingInfo = Info(Score_start:length(Info));
                        AllComaPos = regexp(RemainingInfo,',');
                        ComaPos = AllComaPos(1);
                        Score = RemainingInfo(1:ComaPos-1);
                        % %of matches
                        pct_match_num = (N_matches/LenProt)*100;
                        pct_match = num2str(pct_match_num);
                        t_gw = text(tc_First,Pos+0.05,strcat('GW-TC-Score:',Score,{' '},pct_match,'% id',{' '},'Al-SeCys:',SelCysMatch));
                         
                        if pct_match_num <= 30
                        set(t_gw,'color','r')
                        end 
                        if pct_match_num > 30 & pct_match_num <= 60
                        set(t_gw,'color','m')
                        end 
                        if pct_match_num > 60 & pct_match_num < 90
                        set(t_gw,'color','b')
                        end
                        if pct_match_num >= 90
                        set(t_gw,'color','g')
                        end
                    end
                    end
                end
                end
            end
        end
        
        SumToY = SumToY + 1.5;
        end

        end
    end
    title(Folder)
%         ylim([-2 2])
%          xlim([min(ss_starts) max(ss_ends)])
    xlabel('Contig Position')
    ylabel('Contig')
    saveas(fig,strcat(Folder{1},'_analysis_definitiva'))    
    close;
    end
    end
end
