

In [13]:

    
%matplotlib inline
import nibabel as nib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import metrics
from os.path import join
from nilearn.masking import apply_mask
from nilearn.image import threshold_img,math_img,mean_img
from nilearn import plotting, input_data, image
from scipy.stats import ttest_ind

rsn10 = r'C:\github\a_htywork\20170315_autumnROC\PNAS_Smith09_rsn10.nii.gz'
face_dir = r'C:\expdata\ABIDE_2017_drz5000\FACE_drzstat_N1000'
abide_mask_dir = r'C:\expdata\ABIDE_2017_drz5000\ABIDE12_drz5000_mask'
result_dir = r'C:\temp'



def getAuc_P(drz_ff,mask,target):
    try:
        roivalue = np.mean(apply_mask(drz_ff, mask),1)
        fpr, tpr, thresholds = metrics.roc_curve(target, roivalue)
        auc = metrics.auc(fpr, tpr)
        if auc < 0.5:
            fpr, tpr, thresholds = metrics.roc_curve(target, -roivalue)
            auc = metrics.auc(fpr, tpr)
            
        t, p = ttest_ind(roivalue[(target==0).nonzero()[0]],
                             roivalue[(target==1).nonzero()[0]],
                             equal_var=True)

    except BaseException as e:
        roivalue = target*0        
        auc = 0
        t = 0
        p = 1
        
    return roivalue.reshape(-1,1), auc, t, p



In [18]:

    
face_target = np.ones((43,1))
face_target[:23] = 0
zthr = 4
zstr = str(zthr).replace('.','p')
RSNstr = ['PVN','OPN','LVN','DMN','CEB','MOT','AUD','ECN','RFP','LFP']

'''
with open(join(abide_mask_dir,'subjdx.txt'), 'r') as f:
    abide_target = [int(line.strip()) for line in f]
'''
for ii in range(10):
    y = face_target
    facedata_ff = join(face_dir,'rsn%d' % ii,'allRmaps.nii.gz')
    abide_ff = join(abide_mask_dir,'rsn%d' % ii,'allRmaps.nii.gz')
    
    # get RSN mask with Z threshold
    rsnmask =  math_img('img1>=%f' % zthr,img1=image.index_img(rsn10, ii))
    
        
    abide_maskt1_ff =join(abide_mask_dir,'rsn%d' % ii,'twosample_tfce_corrp_tstat1.nii.gz')
    abide_maskt2_ff =join(abide_mask_dir,'rsn%d' % ii,'twosample_tfce_corrp_tstat2.nii.gz')

    abide_maskt1_ff = math_img('img1>=%f' % 0.95,img1=abide_maskt1_ff)
    abide_maskt2_ff = math_img('img1>=%f' % 0.95,img1=abide_maskt2_ff)
    
    abide_maskt1_ff = math_img('img1*img2',img1=rsnmask, img2=abide_maskt1_ff)  
    abide_maskt2_ff = math_img('img1*img2',img1=rsnmask, img2=abide_maskt2_ff)
    math_img('img1*1.0',img1=abide_maskt1_ff).to_filename(join(result_dir,'A%02dmask.nii.gz' % (ii+1)))
    math_img('img1*1.0',img1=abide_maskt2_ff).to_filename(join(result_dir,'A%02dmask.nii.gz' % (ii+11)))

    
    roivalueR, aucR, tR, pR = getAuc_P(facedata_ff,rsnmask,face_target)
    roivalueA1, aucA1, tA1, pA1 = getAuc_P(facedata_ff,abide_maskt1_ff,face_target)
    roivalueA2, aucA2, tA2, pA2 = getAuc_P(facedata_ff,abide_maskt2_ff,face_target)
    
    plotting.plot_roi(abide_maskt1_ff,title='RSN%d:t1,%d,%.3f' % (ii+1,zthr,pA1),
                      output_file = join(result_dir,'RSN%d_t1_%s.jpg' % (ii+1,zstr)))
    plotting.plot_roi(abide_maskt2_ff,title='RSN%d:t2,%d,%.3f' % (ii+1,zthr,pA2),
                      output_file = join(result_dir,'RSN%d_t2_%s.jpg' % (ii+1,zstr)))



    if ii == 0:
        Rmatrix = roivalueR
        A1matrix = roivalueA1
        A2matrix = roivalueA2
    else:
        Rmatrix = np.hstack((Rmatrix,roivalueR))
        A1matrix = np.hstack((A1matrix,roivalueA1))
        A2matrix = np.hstack((A2matrix,roivalueA2))
    

    #maskt1pxs = abide_maskt1_ff.get_data().sum()
    #maskt2pxs = abide_maskt2_ff.get_data().sum()
    def gpstr(p):
        pstr = '%.3f' % p
        if p < 0.05:pstr += '*'
        if p < 0.01:pstr += '*'
        return pstr
    print('RSN%02d_%s AUC: %.3f(%s)\t t1:%.3f(%s)\t t2:%.3f(%s) ' % (ii+1,RSNstr[ii],aucR,gpstr(pR),aucA1,gpstr(pA1),aucA2,gpstr(pA2)))


np.savetxt(join(result_dir,'tsall_%s.csv' % zstr),np.hstack((face_target,Rmatrix,A1matrix,A2matrix)),fmt='%.5f', delimiter=',')









    



RSN01_PVN AUC: 0.533(0.875)	 t1:0.500(0.931)	 t2:0.593(0.188) 






    



C:\Anaconda3\lib\site-packages\numpy\ma\core.py:4185: UserWarning: Warning: converting a masked element to nan.
  warnings.warn("Warning: converting a masked element to nan.")






    



RSN02_OPN AUC: 0.583(0.637)	 t1:0.000(1.000)	 t2:0.567(0.220) 
RSN03_LVN AUC: 0.593(0.176)	 t1:0.607(0.155)	 t2:0.617(0.204) 
RSN04_DMN AUC: 0.609(0.407)	 t1:0.757(0.010*)	 t2:0.639(0.069) 
RSN05_CEB AUC: 0.567(0.091)	 t1:0.000(1.000)	 t2:0.628(0.036*) 
RSN06_MOT AUC: 0.565(0.313)	 t1:0.604(0.155)	 t2:0.567(0.229) 
RSN07_AUD AUC: 0.628(0.124)	 t1:0.641(0.231)	 t2:0.570(0.304) 
RSN08_ECN AUC: 0.630(0.131)	 t1:0.617(0.135)	 t2:0.663(0.037*) 
RSN09_RFP AUC: 0.622(0.041*)	 t1:0.546(0.685)	 t2:0.663(0.068) 
RSN10_LFP AUC: 0.637(0.160)	 t1:0.667(0.134)	 t2:0.674(0.094)