<jupyter_start><jupyter_text>scRFE Starter CodeV3<jupyter_code># 2020-05-06
# getLabels Function, adding makeOneClassForest Function 
# Imports 
import numpy as np
import pandas as pd
import scanpy as sc
from anndata import read_h5ad
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.feature_selection import SelectFromModel
from sklearn.metrics import accuracy_score
from sklearn.feature_selection import RFE
from sklearn.feature_selection import RFECV
# read in data
adata = read_h5ad('/Users/madelinepark/Downloads/Kidney_facs.h5ad')
sc.pp.normalize_total(adata, target_sum = 1e6, inplace = True)
# adata = read_h5ad('/Users/madelinepark/Downloads/Heart_droplet.h5ad')
adata.shape
# goal: get labels on a per class basis that will go into randomForest function for y
# THE ONE I HAVE BEEN USING THAT WORKS 
def getLabels (dataMatrix, classOfInterest): 
    labelsDict = {}
    for label in np.unique(adata.obs[classOfInterest]):
        lists = []        
        for obs in adata.obs[classOfInterest]:
            if obs == label: 
                lists.append('A')
            else:
                lists.append('B')
        labelsDict[label] = lists #this is usually in line w if and else    
    return labelsDict
# getLabels Unit Test Cases 
# # data with one class, two classes, and two+ classes 
# # droplet vs facs? maybe this doesn't matter 
# CREATE A TEST TO SEE SIMILARITY BETWEEN LISTS (THEY SHOULD NOT BE THE SAME), 
# could use a BOOL?: 
labelsDict = getLabels(dataMatrix = adata, classOfInterest = 'age')
labelsDict['24m']==labelsDict['3m']
# kidneyFacs_getLabels_unitTest = getLabels(dataMatrix = '/Users/madelinepark/Downloads/Kidney_facs.h5ad', classOfInterest = 'age')
# heartDroplet_getLabels_unitTest = getLabels(dataMatrix = '/Users/madelinepark/Downloads/Heart_droplet.h5ad', classOfInterest = 'age')
def makeOneForest (dataMatrix, classOfInterest, labelOfInterest, nEstimators = 10, randomState = 0,  
                  nJobs = -1, oobScore = True, Step = 0.2, Cv =3): 
    
    print('makeOneForest' + labelOfInterest)
    labelsDict = getLabels(dataMatrix, classOfInterest) 

    feat_labels = dataMatrix.var_names #this is equivalent of the genes
    X = dataMatrix.X
    y = labelsDict[labelOfInterest]

    clf = RandomForestClassifier(n_estimators = nEstimators, random_state = randomState, 
                                 n_jobs = nJobs, oob_score = oobScore)
    selector = RFECV(clf, step = Step, cv = Cv)
    
    print('training...')
    clf.fit(X, y)
    selector.fit(X, y)
    feature_selected = feat_labels[selector.support_] 

    return feature_selected, selector.estimator_.feature_importances_ 
FORESTbcell = makeOneForest (dataMatrix = adata, classOfInterest = 'cell_ontology_class',
                             labelOfInterest = 'B cell')
FORESTtcell = makeOneForest (dataMatrix = adata, classOfInterest = 'cell_ontology_class', 
                             labelOfInterest = 'T cell')
FORESTmacro = makeOneForest (dataMatrix = adata, classOfInterest = 'cell_ontology_class',
                                labelOfInterest = 'macrophage')
FORESTmacro[0] == FORESTfencell[0]
good bc the lengths should not be the same
FOREST3 = makeOneForest (dataMatrix = adata, classOfInterest = 'age',
                         labelOfInterest = '3m')
FOREST24 = makeOneForest (dataMatrix = adata, classOfInterest = 'age',
                          labelOfInterest = '24m')
FOREST3[1] == FOREST24[1]
def resultWrite (dataMatrix, classOfInterest, labelOfInterest, nEstimators = 10, randomState = 0,  
                  nJobs = -1, oobScore = True, Step = 0.2, Cv =3):
    print ('result writing')
    results_df = pd.DataFrame() 

    MakeOneForest = makeOneForest(dataMatrix, classOfInterest, labelOfInterest, 
                                  nEstimators, randomState,  nJobs, 
                                  oobScore, Step, Cv)

    feature_selected = MakeOneForest[0]
    feature_gini = MakeOneForest[1]
    
    column_headings = []
    column_headings.append(labelOfInterest)
    column_headings.append(labelOfInterest + '_gini')
    resaux = pd.DataFrame(columns = column_headings)
    resaux[labelOfInterest] = feature_selected
    resaux[labelOfInterest + '_gini'] = (feature_gini)
    results_df = pd.concat([results_df,resaux],axis=1)
    
    return results_df 
def scRFE(dataMatrix, classOfInterest,  nEstimators = 10, randomState = 0,  
                  nJobs = -1, oobScore = True, Step = 0.2, Cv =3):
    
    labelsDict = getLabels(dataMatrix, classOfInterest) 
    
    for label in np.unique(dataMatrix.obs[classOfInterest]):    
        print( 'scRFE' + label)
        
        makeOneForest(dataMatrix = dataMatrix, classOfInterest = classOfInterest, 
                      labelOfInterest = label, nEstimators = nEstimators, 
                      randomState = randomState,  nJobs = nJobs, 
                      oobScore = oobScore, Step = Step, Cv = Cv)
        ResultWrite = resultWrite (dataMatrix = dataMatrix, classOfInterest = classOfInterest, 
                     labelOfInterest = label)
        return ResultWrite
# running with 1000 trees kidney facs data to compare to R 
# scRFE (dataMatrix = adata, classOfInterest = 'cell_ontology_class')
try1 = scRFE (dataMatrix = adata, classOfInterest = 'age')
try1<jupyter_output><empty_output><jupyter_text>Running kidney facs celltype 1000 trees, separate helper functions NOT scRFE to compare with R<jupyter_code># for Sevahn
# 05/13/20
# cells to use: 


# kidney collecting duct pr.3
# kidney collecting duct ep.6

set(adata.obs['cell_ontology_class'])
def getLabels (dataMatrix, classOfInterest): 
    labelsDict = {}
    for label in np.unique(adata.obs[classOfInterest]):
        lists = []        
        for obs in adata.obs[classOfInterest]:
            if obs == label: 
                lists.append('A')
            else:
                lists.append('B')
        labelsDict[label] = lists #this is usually in line w if and else    
    return labelsDict
def makeOneForest_resultWrite (dataMatrix, classOfInterest, labelOfInterest, 
                               nEstimators = 1000, randomState = 0,  nJobs = -1, 
                               oobScore = True, Step = 0.2, Cv =3): 
    
    print('makeOneForest' + labelOfInterest)
    labelsDict = getLabels(dataMatrix, classOfInterest) 

    feat_labels = dataMatrix.var_names #this is equivalent of the genes
    X = dataMatrix.X
    y = labelsDict[labelOfInterest]

    clf = RandomForestClassifier(n_estimators = nEstimators, random_state = randomState, 
                                 n_jobs = nJobs, oob_score = oobScore)
    selector = RFECV(clf, step = Step, cv = Cv)
    
    print('training...')
    clf.fit(X, y)
    selector.fit(X, y)
    feature_selected = feat_labels[selector.support_] 

    print ('result writing')
    feature_gini = selector.estimator_.feature_importances_ 
    
    results_df = pd.DataFrame() 

    column_headings = []
    column_headings.append(labelOfInterest)
    column_headings.append(labelOfInterest + '_gini')
    resaux = pd.DataFrame(columns = column_headings)
    resaux[labelOfInterest] = feature_selected
    resaux[labelOfInterest + '_gini'] = (feature_gini)
    results_df = pd.concat([results_df,resaux],axis=1)
    
    return results_df
kidney1000TCELL = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'T cell')
kidney1000macrophage = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'macrophage')
kidney1000fenestrated = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'fenestrated cell')
kidney1000epithelialProximal = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'epithelial cell of proximal tubule')
kidney1000loopHenle = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'kidney loop of Henle ascending limb epithelial cell')
kidney1000collectingDuctEP  = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'kidney collecting duct epithelial cell')
kidney1000collectingDuctPrincipal  = makeOneForest_resultWrite(dataMatrix = adata, 
                                            classOfInterest = 'cell_ontology_class', 
              labelOfInterest = 'kidney collecting duct principal cell')
cd desktop
# save all as csv to desktop before I lose them! 
# kidney1000collectingDuctPrincipal.to_csv('kidney1000collectingDuctPrincipal.csv',index=False)
# kidney1000collectingDuctEP.to_csv('kidney1000collectingDuctEP.csv',index=False)
# kidney1000loopHenle.to_csv('kidney1000loopHenle.csv', index = False)
# kidney1000epithelialProximal.to_csv('kidney1000epithelialProximal.csv', index = False)
# kidney1000fenestrated.to_csv('kidney1000fenestrated.csv', index = False)
# kidney1000macrophage.to_csv('kidney1000macrophage.csv', index = False)
# kidney1000TCELL.to_csv('kidney1000TCELL.csv', index = False)
# newkidney = kidney1000TCELL.sort_values(by= 'T cell_gini', ascending = False)
five = kidney1000macrophage.sort_values(
by = 'macrophage_gini', ascending = False)
five.to_csv('kidney1000macrophageSORTED.csv', index = False)
set(adata.obs['cell_ontology_class'])
# is it printing same results bc of how I am looping? like it's getting over ridden?
def makeOneClassForestOld(dataMatrix, classOfInterest, step=0.2, cv=3,  n_estimators=1000, random_state=0, n_jobs=-1, oob_score=True, min_samples_split=2, min_samples_leaf=1):
#     results_df = pd.DataFrame() 
    results_df = pd.DataFrame() 
    labelsDict = getLabels(dataMatrix, classOfInterest) 
    print(labelsDict.keys())
    for label in list(labelsDict.keys())[0:3]:
        print(label)
        feat_labels = dataMatrix.var_names #this is equivalent of the genes
#         X = dataMatrix.X
#         y1 = np.array([labelsDict[label]]) #series of strings, converted from list
# #         y = y1.transpose() #maybe its not transpose, but RESHAPE
#         y = y1.reshape(1400,1)
        X = dataMatrix.X
        y = labelsDict[label]
#         print(type(X))
#         print((y.shape))
        clf = RandomForestClassifier(n_estimators, random_state, n_jobs, oob_score, min_samples_split, min_samples_leaf)
        selector = RFECV(clf, step, cv) #when I hardcode args for this and clf, it does not error
        
        clf.fit(X, y)
        selector.fit(X, y)
        feature_selected = feat_labels[selector.support_] 
        
        print('result writing')
#         column_headings = []
#         column_headings.append(c)
#         column_headings.append(c + '_gini')
#         resaux = pd.DataFrame(columns=column_headings)
#         resaux[c] = feature_selected
#         resaux[c + '_gini'] = (selector.estimator_.feature_importances_)

#         print(feature_selected)
#         print (selector.estimator_.feature_importances_)
#         results_df = pd.concat([results_df,resaux],axis=1)
#     results_df 
# call getLabels
labelsDict = getLabels(dataMatrix = adata, classOfInterest = 'age')

values = labelsDict.values()
values_list = list(values)
len(values_list[0]) #so this works 
len(labelsDict['24m'])
dataMatrix = adata
classOfInterest = 'age'
labelsDict = getLabels(dataMatrix, classOfInterest) 
results_df = pd.DataFrame() 

for label in list(labelsDict.keys()):
    print(label)
    feat_labels = dataMatrix.var_names #this is equivalent of the genes
    X = dataMatrix.X
    y = labelsDict[label] #dim = (1,1400)
#     print(len(y))
#     print(type(y))
#     y = np.array([labelsDict[label]]) #dim = (1,1400)
#     y.shape = (1400,1)
#     print(y.shape)
#     print(type(y))
#     print(type(X))
#     print((y.shape))

#     #problem is that the length of y1 is 1 but the shape is 1400,1
    print('training...')
    clf = RandomForestClassifier(n_estimators=10, random_state=0, n_jobs=-1, oob_score=True)
    selector = RFECV(clf, step=0.2, cv=3) # step = % rounded down at each iteration
    clf.fit(adata.X, y) #confused if this y is labels or the genes because feat_labels is genes
    selector.fit(X, y)
    feature_selected = feat_labels[selector.support_] 

#     print('result writing...')
#     column_headings = []
#     column_headings.append(c)
#     column_headings.append(c + '_gini')
#     resaux = pd.DataFrame(columns=column_headings)
#     resaux[c] = feature_selected
#     resaux[c + '_gini'] = (selector.estimator_.feature_importances_)

    print(feature_selected)
    print (selector.estimator_.feature_importances_)
#     results_df = pd.concat([results_df,resaux],axis=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09446179629770701, 'fraction_numerical': 0.018054391711739164, 'mean_word_length': 3.097378277153558, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '216509', 'license': 'MIT', 'n_tokens_mistral': 4374, 'n_tokens_neox': 3989, 'n_words': 917, 'path': 'scripts/practiceScripts/scRFE_starter_codeV3.ipynb', 'repo_name': 'czbiohub/scRFE', 'size': '13127'}