function [mNet,mNet_scrub,mNet_scrub_GSR,IOR] = grabMRI(TC_fullfile,GSR,scratch,rate_data)


[matrix] = load(TC_fullfile);
    funcmat         = double(matrix.matrix);
    ts.ts           = funcmat;
    ts.tr           = 2.48;   
    ts.Nnodes       = size(funcmat,2);
    ts.Ntimepoints  = size(funcmat,1);
    ts.Nsubjects    = 1;
    ts.NnodesOrig   = ts.Nnodes;
    ts.NtimepointsPerSubject = size(funcmat,1);
    ts.DD           = 1:ts.Nnodes;
    ts.UNK          = [];   

netmats       = nets_netmats_md(ts,1,0,'corr'); % correlation transformed 2 Z-scores
mNet          = reshape(netmats(1,:),[sqrt(length(netmats)) sqrt(length(netmats))]);
[IOR_mNet] = InOutRatio(mNet);
% Scrubbing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ts_scrub = ts;
    if ~isempty(scratch)
        funcmat_scrub = ts_scrub.ts;
        funcmat_scrub(scratch,:) = [];
        ts_scrub.ts = funcmat_scrub;
        ts_scrub.Ntimepoints  = size(funcmat_scrub,1);
        ts_scrub.NtimepointsPerSubject = size(funcmat_scrub,1);
    else
        funcmat_scrub = ts_scrub.ts;
    end
netmats_scrub       = nets_netmats_md(ts_scrub,1,0,'corr'); % correlation transformed 2 Z-scores
mNet_scrub          = reshape(netmats_scrub(1,:),[sqrt(length(netmats_scrub)) sqrt(length(netmats_scrub))]);
IOR_mNet_scrub = InOutRatio(mNet_scrub); 
% Global signal regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if exist("GSR","var")  
     if ~isempty(scratch)
     GSR(scratch)=[];
     end
     yadj = zeros(size(ts_scrub.ts));
     X = detrend(GSR,'constant');
     for v = 1:ts.Nnodes
         [Q, R]=qr(X,0);
         b = R\(Q'*funcmat_scrub(:,v));
         yvoxmodel = b(1)*X(:,1);

    % Adjust time-series
         yadj(:,v) = funcmat_scrub(:,v) - yvoxmodel;
     end
    ts_scrub.ts = double(yadj);
 end
netmats_scrub_GSR       = nets_netmats_md(ts_scrub,1,0,'corr'); % correlation transformed 2 Z-scores
mNet_scrub_GSR          = reshape(netmats_scrub_GSR(1,:),[sqrt(length(netmats_scrub_GSR)) sqrt(length(netmats_scrub_GSR))]);
IOR_mNet_scrub_GSR = InOutRatio(mNet_scrub_GSR);

IOR = IOR_mNet_scrub_GSR;

if rate_data == 1
 % compute matrix plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure('units','normalized','outerposition',[0 0 1 1])
         
    subplot(131)
        plotmat_values_rsfMRI(mNet,[],[-1 1],0,[]);
        hold on
        axis square;
            plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
            plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title(strcat('raw mNet (IOR: ',num2str(round(IOR_mNet,2)),')'))
    subplot(132)
        plotmat_values_rsfMRI(mNet_scrub,[],[-1 1],0,[]);
        hold on
        axis square;
            plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
            plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title({strcat('mNet scratched (IOR: ',num2str(round(IOR_mNet_scrub,2)),')')})
    subplot(133)
        plotmat_values_rsfMRI(mNet_scrub_GSR,[],[-1 1],0,[]);
        hold on
        axis square;
            plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
            plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title(strcat('mNet scratched + GSR (IOR: ',num2str(round(IOR_mNet_scrub_GSR,2)),')'))
end
end