v2.15

Author: PhilippSpahn

Scripts/AlignReads.py

@@ -360,6 +360,9 @@ def MapAndCount(sample):
     for k in range(L):
         GuideCounts.write(str(sgIDs[k]) + '\t'+ str(geneIDs[k]) + '\t' + str(ReadsPerGuide[k]) + '\n')
     GuideCounts.close()
+    # No-mapping warning
+    if sum(ReadsPerGuide) == 0:
+        print('!! ERROR: Zero total read counts! Check library file and index !!')
     # Read counts per gene in library       
     print('Counting reads per gene ...')   
     global GeneList

Scripts/FindHits.py

@@ -157,10 +157,10 @@ def PrepareHitList(sample):
                                      'control stdev [norm.]': [numpy.sqrt(sigma2[k]) for k in range(L)],
                                      'fold change': [fc[k] for k in range(L)],   
                                      'p-value': ['%.2E' % Decimal(NBpval[k]) for k in range(L)],
-                                     'adj. p-value': ['%.2E' % Decimal(NBpval_0[k]) for k in range(L)],                                                 
+                                     'FDR': ['%.2E' % Decimal(NBpval_0[k]) for k in range(L)],                                                 
                                      'significant': [str(significant[k]) for k in range(L)]},
                             columns = ['sgRNA','gene','counts [norm.]','control mean [norm.]',\
-                            'control stdev [norm.]','fold change','p-value','adj. p-value','significant'])
+                            'control stdev [norm.]','fold change','p-value','FDR','significant'])
     if ScreenType == 'enrichment':
         Results_df_0 = Results_df.sort_values(['significant','fold change'],ascending=[0,0])                
     elif ScreenType == 'depletion':

Scripts/PinAPL.py

@@ -170,8 +170,10 @@
     StatMsg = 'Time elapsed [hours]: ' + '%.3f' % time_elapsed +'\n'
     os.system('python -u PrintStatus.py TimeStamp "'+StatMsg+'" 2>&1 | tee -a PinAPL-Py.log')
 
-# Move Log File 
+# Move Log Files 
 if not os.path.exists(LogFileDir):
 	os.makedirs(LogFileDir) 
 os.system('cp configuration.yaml '+LogFileDir)
-os.system('cp PinAPL-Py.log '+LogFileDir)
+os.system('cp PinAPL-Py.log '+LogFileDir)
+os.chdir(WorkingDir)
+os.system('cp DataSheet.xlsx '+LogFileDir)

Scripts/PlotCounts.py

@@ -120,9 +120,9 @@ def GOI_Scatterplot(sample,GOI='None'):
         if svg:
             plt.savefig(sample+' '+GOI+' counts.svg')
     else:
-        plt.savefig(sample+' '+' counts.png', dpi=res)
+        plt.savefig(sample+' counts.png', dpi=res)
         if svg:
-            plt.savefig(sample+' '+' counts.svg')
+            plt.savefig(sample+' counts.svg')
     plt.close()
 
     # ------------------------------------------------

Scripts/RankGenes.py

@@ -85,7 +85,7 @@ def compute_aRRA(g):
     GOI = geneList[g]    
     GOI_ranks_sig = list()    
     for i in range(L):
-        if GOI == genes[i] and (NB_pval[i] < alpha):
+        if GOI == genes[i] and (NB_pval[i] < P_0):
             GOI_ranks_sig.append(ranks[i])
     j = len(GOI_ranks_sig)
     if j>0:
@@ -103,7 +103,7 @@ def compute_aRRA(g):
 def compute_aRRA_null(I):
     I_ranks_sig = list()
     for i in range(L):
-        if (i in I) and (NB_pval[i] < alpha):
+        if (i in I) and (NB_pval[i] < P_0):
             I_ranks_sig.append(ranks[i])
     j = len(I_ranks_sig)
     if j>0:
@@ -152,7 +152,7 @@ def GeneRankingAnalysis(sample):
     ListDir = config['HitDir']
     EffDir = config['EffDir']
     GeneDir = config['GeneDir']
-    global alpha; alpha = config['alpha']            
+    alpha = config['alpha']            
     padj = config['padj']
     num_cores = multiprocessing.cpu_count()
     global r; r = config['sgRNAsPerGene']
@@ -162,6 +162,7 @@ def GeneRankingAnalysis(sample):
     pvalDir = config['pvalDir']
     res = config['dpi']
     svg = config['svg']
+    global P_0; P_0 = 0.05 # p-value Threshold for aRRA analysis
 
     # ------------------------------------------------
     # Read sgRNA enrichment/depletion table
@@ -224,7 +225,7 @@ def GeneRankingAnalysis(sample):
             plt.savefig(sample+'_sgRNA_Efficiency.svg')               
     else: # no control replicates
         print('WARNING: No control replicates found! No significant sgRNAs counted.')
-        sigGuides = ['N/A' for k in range(len(sigGuides))]
+        sigGuides = ['N/A' for k in range(G)]
         # Plot histogram
         if not os.path.exists(EffDir):
             os.makedirs(EffDir)
@@ -316,6 +317,7 @@ def GeneRankingAnalysis(sample):
             SortFlag = True
             metric = [-1 for k in range(G)]
             metric_pval = [-1 for k in range(G)]
+            metric_pval0 = [-1 for k in range(G)]            
             metric_sig = ['N/A' for k in range(G)]
     elif GeneMetric == 'STARS':
         # -------------------------------------------------        
@@ -376,8 +378,9 @@ def GeneRankingAnalysis(sample):
     # -------------------------------------------------  
     # p-value plots
     # -------------------------------------------------  
-    print('Plotting p-values ...')
-    pvalHist_metric(metric_pval,metric_pval0,GeneMetric,pvalDir,sample,res,svg)           
+    if min(NB_pval) > -1:    
+        print('Plotting p-values ...')
+        pvalHist_metric(metric_pval,metric_pval0,GeneMetric,pvalDir,sample,res,svg)           
 
     # -------------------------------------------------  
     # Output list
@@ -389,15 +392,16 @@ def GeneRankingAnalysis(sample):
     Results_df = pandas.DataFrame(data = {'gene': [geneList[g] for g in range(G)],
                                     GeneMetric: [metric[g] for g in range(G)],
                                      GeneMetric+' p_value': ['%.2E' % Decimal(metric_pval[g]) for g in range(G)],
-                                    GeneMetric+' adj. p_value': ['%.2E' % Decimal(metric_pval0[g]) for g in range(G)],
-                                     'significant': [metric_sig[g] for g in range(G)],                                                    
+                                    GeneMetric+' FDR': ['%.2E' % Decimal(metric_pval0[g]) for g in range(G)],
+                                     'significant': [str(metric_sig[g]) for g in range(G)],                                                    
                                      '# signif. sgRNAs': [sigGuides[g] for g in range(G)]},
-                            columns = ['gene',GeneMetric,GeneMetric+' p_value',GeneMetric+' adj. p_value',\
+                            columns = ['gene',GeneMetric,GeneMetric+' p_value',GeneMetric+' FDR',\
                             'significant','# signif. sgRNAs'])
     Results_df_0 = Results_df.sort_values(['significant',GeneMetric],ascending=[False,SortFlag])
     GeneListFilename = filename[0:-14]+'_'+GeneMetric+'_'+'P'+str(Np)+'_GeneList.tsv'
     Results_df_0.to_csv(GeneListFilename, sep = '\t', index = False)      
     if SheetFormat == 'xlsx':
+        print('Converting to xlsx ...')
         GeneListFilename = filename[0:-14]+'_'+GeneMetric+'_'+'P'+str(Np)+'_GeneList.xlsx'
         Results_df_0.to_excel(GeneListFilename)              
 

configuration.yaml

@@ -9,7 +9,7 @@
 ScreenType: 'enrichment'                    # type of screen ['enrichment'/'depletion']
 
 # LIBRARY PARAMETERS
-LibFilename: 'GeCKOv2_library.tsv'            # filename of library spreadsheet  
+LibFilename: 'GeCKOv2_Human.tsv'            # filename of library spreadsheet  
 seq_5_end: 'TCTTGTGGAAAGGACGAAACACCG'       # sequence 5' of sgRNA in read
 seq_3_end: 'GTTTTAGAGCTAGAAATAGCAAGTT'       # sequence 3' of sgRNA in read
 NonTargetPrefix: 'NonTargeting'        # prefix for non-targeting sgRNAs in library (keep at default if none present)
@@ -22,9 +22,9 @@ keepCutReads: False                    # keep files containing trimmed reads? ['
 Normalization: 'cpm'                  # Normalization of read counts (counts per mio reads / size-factor) ['cpm'/'size']
 Cutoff: 0                           # set counts (cpm) below this number to 0 
 GeneMetric: 'aRRA'                    # metric for gene ranking ['aRRA'/'STARS'/'ES']
-ClusterBy: 'variance'                  # clustering criterion ['variance'/'counts']
+ClusterBy: 'counts'                  # clustering criterion ['variance'/'counts']
 HitListFormat: 'tsv'                  # Format of results spreadsheets (sgRNA hits and gene ranking) ['tsv'/'xlsx']
-alpha: 0.01                            # significance level    
+alpha: 0.10                            # false discovery rate
 padj: 'fdr_bh'                        # method for p-value adjustment ['fdr_bh'/'fdr_tsbh']