timydaley / CRISPhieRmix

8 stars 5 forks source link

I am getting an error that halts execution of the script. #1

Closed sbodapati closed 5 years ago

sbodapati commented 5 years ago

Error is "Error in if (abs(loglike - prevloglike)/n_obs < tol | iter > max_iter) { : missing value where TRUE/FALSE needed"

The input is CRISPhieRmix_Simulation_l2fc.txt control is Simulation_l2fc_control.txt

CRISPhieRmix_Simulation_l2fc.txt

timydaley commented 5 years ago

Can you give the exact commands that produce the error so that I can recreate the error?

sbodapati commented 5 years ago

f_name = 'Simulation' dataset = 'Simulation' controlGuideName = 'Simulation_l2fc_control.csv'

l2f_filename = paste("/CRISPRBenchmarking/input/CRISPhieRmix",f_name,"_l2fc",".csv", sep='') l2fc = read.table(file = l2f_filename, header = TRUE, sep=',') log2fc = l2fc$x print(length(log2fc))

filename = paste("CRISPRBenchmarking/input/CRISPhieRmix",f_name,".csv", sep='') testCounts = read.table(file = filename, header = TRUE, sep=',') geneID = testCounts$gene_id geneID = factor(geneID,levels=unique(geneID)) print(length(geneID))

controlGuides_filename = paste("CRISPR_Benchmarking/controls/",controlGuideName, sep='') controlGuides = read.table(file = controlGuides_filename, header = TRUE, sep=',') control2fc = controlGuides$x

CRISPhieRmixResults = CRISPhieRmix::CRISPhieRmix(x = log2fc, geneIds = geneID, negCtrl = control2fc, BIMODAL=TRUE)

timydaley commented 5 years ago

The geneIDs were not attached. Please generate a single file containing the log2fc and the geneIDs as well as a minimal script to reproduce the error.

sbodapati commented 5 years ago

Here are all of the input files and the script to run them.

[CRISPhieRmix_Simulation.txt] = first col is the gene names that correspond to the gene targeting guide l2fc below. (https://github.com/timydaley/CRISPhieRmix/files/3086292/CRISPhieRmix_Simulation.txt) [CRISPhieRmix_Simulation_l2fc.txt] = gene targeting l2fc (https://github.com/timydaley/CRISPhieRmix/files/3086293/CRISPhieRmix_Simulation_l2fc.txt) [Simulation_l2fc_control.txt] = Control l2fc (https://github.com/timydaley/CRISPhieRmix/files/3086305/Simulation_l2fc_control.txt) [CRISPhieRmix_test.txt] = R Script (https://github.com/timydaley/CRISPhieRmix/files/3086309/CRISPhieRmix_test.txt)

timydaley commented 5 years ago

Adding the following to the script results in no failure: if(!exists("CRISPhieRmixResults")){

failure CRISPhieRmix, run in unimodal mode

if (length(control2fc)>0){ CRISPhieRmixResults = CRISPhieRmix::CRISPhieRmix(x = log2fc, geneIds = geneID, negCtrl = control2fc, mu = -4, screenType = "LOF") # geneRanks = CRISPhieRmixResults$FDR } else { CRISPhieRmixResults = CRISPhieRmix::CRISPhieRmix(x = log2fc, geneIds = geneID, mu=-4) # geneRanks = CRISPhieRmixResults$FDR } }

}

The issue is that there really is no signal here. Attached is a figure of gene targetting vs control densities, and the control is wider than the gene targetting, indicating that there is no signal. I would expect any reasonable method to fail to find any significant genes with this data.
test_data_control_vs_gene.pdf

Basically, I'm comfortable with the behavior of CRISPhieRmix with this data.