I'm running into an error on the MAGeCK-FLUTE step, it seems to be complaining about genes that have a duplicate Entrez ID? My command was

nextflow run nf-core/crisprseq -r dev -c /pipelines/nfcore/nextflow.config -work-dir /nextflow/work --input crisprseq_samplesheet.csv --outdir nextflow_out --analysis screening --library gRNA_annotated_geneSymbol_withID.tsv --contrasts crisprseq_contrasts.csv --reference_fasta GRCh38.primary_assembly.genome.fa

contrasts are set up as

reference;treatment
control;treated_5hr
control;treated_7hr

any clues?

ERROR ~ Error executing process > 'NFCORE_CRISPRSEQ:CRISPRSEQ_SCREENING:MAGECK_FLUTEMLE_CONTRASTS (treated_5hr_vs_control)'

Caused by:
  Process `NFCORE_CRISPRSEQ:CRISPRSEQ_SCREENING:MAGECK_FLUTEMLE_CONTRASTS (treated_5hr_vs_control)` terminated with an error exit status (1)

Command executed [/homes8/jsimon/.nextflow/assets/nf-core/crisprseq/templates/template_fluteMLE.R]:

  #!/usr/bin/env Rscript
  #### author: Laurence Kuhlburger
  #### Released under the MIT license. See git repository (<https://github.com/nf-core/crisprseq>) for full license text.
  ####
  #### graphs mageck MLE

  library(MAGeCKFlute)
  library(clusterProfiler)
  library(ggplot2)

  library(pathview)
  options(ggrepel.max.overlaps = Inf)
  mle <- read.table("treated_5hr_vs_control.gene_summary.txt", header = TRUE, sep = "   ", stringsAsFactors = FALSE)

  if("treated_5hr_vs_control" == "day0") {
      beta_strings <- grep("\\.beta", colnames(mle), value = TRUE)
      before_beta <- sub("\\.beta.*", "", beta_strings)
      unique_strings <- unique(before_beta)
      for(i in unique_strings) {
          FluteMLE(mle, treatname= i, proj=i, <http://pathview.top|pathview.top>=5)
          }
  } else {
      beta_strings <- grep("\\.beta", colnames(mle), value = TRUE)
      before_beta <- sub("\\.beta.*", "", beta_strings)
      unique_strings <- unique(before_beta)
      for(i in unique_strings) {
          FluteMLE(mle, treatname= i, proj=i,  , <http://pathview.top|pathview.top>=5)
      }
  }

  version_file_path <- "versions.yml"
  version_flute <- paste(unlist(packageVersion("MAGeCKFlute")), collapse = ".")
  version_ggplot <- paste(unlist(packageVersion("ggplot2")), collapse = ".")
  version_clusterprofiler <- paste(unlist(packageVersion("clusterProfiler")), collapse = ".")
  version_pathview <- paste(unlist(packageVersion("pathview")), collapse = ".")

  f <- file(version_file_path, "w")
  writeLines('"NFCORE_CRISPRSEQ:CRISPRSEQ_SCREENING:MAGECK_FLUTEMLE_CONTRASTS":', f)
  writeLines("    MAGeCKFlute: ", f, sep = "")
  writeLines(version_flute, f)
  writeLines("    ggplot2: ", f, sep = "")
  writeLines(version_ggplot, f)
  writeLines("    clusterProfiler: ", f, sep = "")
  writeLines(version_clusterprofiler, f)
  writeLines("    pathview: ", f, sep = "")
  writeLines(version_pathview, f)
  close(f)

Command exit status:
  1

Command output:
  (empty)

Command error:
  36 genes have duplicate Entrez IDs: HIST1H2BC, HIST1H2BG, HIST1H2BF, PLEKHG7, DEC1, HIST1H2BI, MUM1, KIAA1107, AKAP2, NEBL, HN1L, SLC35E2, WWTR1, NOV, DCDC1, C10orf12, ZASP, HYAL1, MICALCL, QARS, LOXL2, CC2D2B, GIF, RAB34, C16orf47, TMEM257, C17orf47, ARHGAP42, MARS, TXNRD3, RRM2, HIST1H2BE, SEPT6, SEM1, C9orf47, DUPD1
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  loading from cache
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  loading from cache
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  loading from cache
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  loading from cache
  2024-07-03 19:14:43 # Enrichment analysis of Positive and Negative selected genes ...
  2024-07-03 19:14:43 # Running KEGG+REACTOME+GOBP+Complex enrichment analysis
        446 genes are mapped ...
  2024-07-03 19:14:52 # Running KEGG+REACTOME+GOBP+Complex enrichment analysis
        330 genes are mapped ...
  2024-07-03 19:15:07 # Starting plot kegg pathways for  Group A
  Info: Downloading xml files for hsa04115, 1/1 pathways..
  Info: Downloading png files for hsa04115, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /jsimonlab/scratch/jsimon/nextflow/work/20/2cd6d926916e086179100326a47cbe
  Info: Writing image file hsa04115.pathview.multi.png
  Info: Downloading xml files for hsa04215, 1/1 pathways..
  Info: Downloading png files for hsa04215, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /jsimonlab/scratch/jsimon/nextflow/work/20/2cd6d926916e086179100326a47cbe
  Info: Writing image file hsa04215.pathview.multi.png
  Info: some node width is different from others, and hence adjusted!
  pathview failed for hsa04215
  Info: Downloading xml files for hsa05222, 1/1 pathways..
  Info: Downloading png files for hsa05222, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /jsimonlab/scratch/jsimon/nextflow/work/20/2cd6d926916e086179100326a47cbe
  Info: Writing image file hsa05222.pathview.multi.png
  Info: Downloading xml files for hsa04210, 1/1 pathways..
  Info: Downloading png files for hsa04210, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /jsimonlab/scratch/jsimon/nextflow/work/20/2cd6d926916e086179100326a47cbe
  Info: Writing image file hsa04210.pathview.multi.png
  Info: Downloading xml files for hsa01524, 1/1 pathways..
  Info: Downloading png files for hsa01524, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /jsimonlab/scratch/jsimon/nextflow/work/20/2cd6d926916e086179100326a47cbe
  Info: Writing image file hsa01524.pathview.multi.png
  Error in png::readPNG(figure, native = FALSE) : file is not in PNG format
  Calls: FluteMLE ... arrangePathview -> lapply -> FUN -> <Anonymous> -> <Anonymous>
  In addition: There were 12 warnings (use warnings() to see them)
  Execution halted

hi Jeremy! Could you confirm that mapping worked with

--reference_fasta GRCh38.primary_assembly.genome.fa

like how does your count table look like from mageck?

both my counts and results from mageck mle look good! i can share snippets if helpful

No shouldn't be the case. So there aren't duplicates in the count table or the mle output one ? Sorry you're the first user that used the reference fasta so I'm also keen to hear the bowtie2 ran for you?

duplicate gene names there, you mean?

Yes :)

just grepped for one of the "duplicated" gene IDs in my mageck count table, looks like this:

sgRNA16220	HIST1H2BC	336	86	102	120	14	105	171	129	72
sgRNA16221	HIST1H2BC	436	366	172	262	330	395	695	420	361
sgRNA16222	HIST1H2BC	346	780	247	392	578	482	431	519	177

Ok i think you can run the command without

--reference_fasta GRCh38.primary_assembly.genome.fa

and it should work the same then

sure thing-

also just a general question, would the mageckflute output be of use to you as a user? I was thinking about making it not run automatically since it seems to be causing issues

if i understand what it does correctly, it sounds potentially useful. my goal here is to compare our treated_5hr vs control, treated_7hr vs control to get genes and pathways significantly enriched in each separately (and in common)

no shouldn’t change anything, here it’d show you enriched pathways against all cell lines in depmap for each comparison

depmap data used for understanding essentiality here?

Yes

yes, human iPS

okay i’ll also have a look if it’s in the depmap database so that you use mageckflute only specifying that particular cell line. I’d then introduce a new param for users along those lines

is it not helpful to test against all cell lines in depmap?

I think it’d give you a broader overview of how your genes/experiments perform against all cell lines? But if you have a treatment and the depmap is only essential genes with no ttt it’d give you a more up close view if you pick one specific cell line?

k - in this case im not 100% sure there are any stem or induced pluripotent cells in depmap but still looking through the list

Yes also don’t think i’ve come accross any

ill have to read up more on how bagel2 works

hm good point. Would you mind opening an issue for that so i keep in mind to improve the doc? I’ll play around with your count table and the day0 label to see what can be obtained, but i think

--day0_label 'Untreated_0hr_Rep1,Untreated_0hr_Rep2,Untreated_0hr_Rep3'

would be the correct way to go for your experiment, haven’t tried it though

let me know if youd like an issue for tracking purposes on the original error posted above ^^

Hi Jeremy, I am trying to reproduce your error with the following command :

nextflow run . -c user.config --count_table ../jeremy/count_table.count.txt --outdir test_jeremy --analysis screening --mle_design_matrix ../jeremy/treated_7hr_vs_control.txt -profile docker -resume

and get another error :

Command error:
  retrieving 1 resource
  loading from cache
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  downloading 1 resources
  retrieving 1 resource
  loading from cache
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  loading from cache
  snapshotDate(): 2022-10-31
  see ?depmap and browseVignettes('depmap') for documentation
  loading from cache
  2024-07-10 12:02:37 # Enrichment analysis of Positive and Negative selected genes ...
  2024-07-10 12:02:37 # Running KEGG+REACTOME+GOBP+Complex enrichment analysis
        484 genes are mapped ...
  2024-07-10 12:02:41 # Running KEGG+REACTOME+GOBP+Complex enrichment analysis
        301 genes are mapped ...
  2024-07-10 12:02:49 # Starting plot kegg pathways for  Group A
  Info: Downloading xml files for hsa04115, 1/1 pathways..
  Info: Downloading png files for hsa04115, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /Users/laurencekuhlburger/Documents/crisprseq/crisprseq/work/d7/ef8f81048ef6d00888210cdd1e6731
  Info: Writing image file hsa04115.pathview.multi.png
  Info: Downloading xml files for hsa04215, 1/1 pathways..
  Info: Downloading png files for hsa04215, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /Users/laurencekuhlburger/Documents/crisprseq/crisprseq/work/d7/ef8f81048ef6d00888210cdd1e6731
  Info: Writing image file hsa04215.pathview.multi.png
  Info: some node width is different from others, and hence adjusted!
  pathview failed for hsa04215
  Info: Downloading xml files for hsa01524, 1/1 pathways..
  Info: Downloading png files for hsa01524, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /Users/laurencekuhlburger/Documents/crisprseq/crisprseq/work/d7/ef8f81048ef6d00888210cdd1e6731
  Info: Writing image file hsa01524.pathview.multi.png
  Info: Downloading xml files for hsa04210, 1/1 pathways..
  Info: Downloading png files for hsa04210, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /Users/laurencekuhlburger/Documents/crisprseq/crisprseq/work/d7/ef8f81048ef6d00888210cdd1e6731
  Info: Writing image file hsa04210.pathview.multi.png
  Info: Downloading xml files for hsa03040, 1/1 pathways..
  Info: Downloading png files for hsa03040, 1/1 pathways..
  'select()' returned 1:1 mapping between keys and columns
  Info: Working in directory /Users/laurencekuhlburger/Documents/crisprseq/crisprseq/work/d7/ef8f81048ef6d00888210cdd1e6731
  Info: Writing image file hsa03040.pathview.multi.png
  Error in png::readPNG(figure, native = FALSE) : file is not in PNG format
  Calls: FluteMLE ... arrangePathview -> lapply -> FUN -> <Anonymous> -> <Anonymous>
  In addition: There were 12 warnings (use warnings() to see them)
  Execution halted

which i think i’ve seen before on pathway enrichment

hm my current setup is via singularity on our hpc. im not 100% sure if/how i can switch

No i don’t think so

Thanks! So I can keep my singularity profile but instead specify -r dev, correct?

yes and also specify

laurencekuhl/crisprseq -r dev

instead of the nf-core one

Will give that a go, thanks!

yay 🙂 thanks a lot for reaching out and opening issues, appreciate it! if you see anything else don’t hesitate to reach out

thanks!

No since you defined contrasts you don’t have to, i never used day0_label either

one last question, i hope! my

mageck count

step shows my sample labels as

--sample-label treated_5hr,treated_5hr,treated_5hr,treated_7hr,treated_7hr,treated_7hr,control,control,control

not as the sample labels (first column) in my sample sheet

Etoposide_5hr_Rep1,Etoposide_5hr_Rep2,Etoposide_5hr_Rep3,Etoposide_7hr_Rep1,Etoposide_7hr_Rep2,Etoposide_7hr_Rep3,Untreated_0hr_Rep1,Untreated_0hr_Rep2,Untreated_0hr_Rep3

and as a result, i can't tell my replicates apart in my counts tables (

count_table.count.txt

,

count_table.count_normalized.txt

) is there an additional metadata column i can pass in my samplesheet, or is there another parameter somewhere ive missed?

sorry Jeremy didn’t have time to have a closer look, hope to get back to you by tomorrow”

no worries - i opened it as an issue on github too in case others are facing a similar issue

I think if you put the sample labels as the condition in your samplesheet as well you’d get what you need eventhough it’s not the best i agree

oh really? i feel like "condition" makes sense in the case of replicates

how about outputting one count table with the sample names and one with the condition names then?

Is my case a rare one? Do most people not have replicates?

no i think it’s quite common

If most are set up like mine then my opinion is the table should be by Sample name, and leave it to the user to know how that maps to conditions

yes because the design matrix needs to have the sample names

Is that how it's working right now? I feel like the design matrix it created for me has only conditions as row names

which right now works because your count table also has conditions as column names but doesn’t allow you to see which sample is which

I think my vote would be to mimic the rnaseq setup here

i’m just not sure if i can get it to output a count table with column names for sample name and then translate it to conditions

Which tool does the gene essentiality test?

mageck mle does gene essentiality

my interpretation of this is in the MAGeCK-MLE docs is that the design matrix rows would be samples, not conditions. is that incorrect?

Yes i think this is an implementation mistake on my end which was why i was thinking of kicking out conditions, it’s just confusing users 😞

how would you change the contrasts specification to match then, if no condition column is supplied?

Hi Jeremy, just added

This label is not mandatory as in case you are running time series. If you wish to run MAGeCK MLE with the day0 label you can do so by specifying `--day0_label` and the sample names that should be used as day0. The contrast will then be automatically adjusted for the other days.

to the day0 label

thanks, that is a helpful add

I have a solution in mind i still need to test out but it might make the pipeline slower 😕 i would run mageck count twice, once with the sample labels and once with the condition labels

how could they be added as labels?

yes they should be, but if you want to be on the safe side i’d suggest rerunning the pipeline and putting the names of your samples in the condition column

okay- do i adjust the contrasts file too then? to something like

reference;treatment
control_1,control_2,control_3;treated_5hr_1,treated_5hr_2,treated_5hr_3
control_1,control_2,control_3;treated_7hr_1,treated_7hr_2,treated_7hr_3

Thanks for reopening, I think since this one is a recurring ask I'll add it as a variable, can I contact you with a dev version as soon as I'm done?

Yes definitely!

could you try with you nextflow run with

-r *control-sgrna*

instead of dev and add

--mle_control_sgrna 'your_file.txt'

?

from your repo, yes?

No it's on nf-core, got a bit tangled up

hmm it seems it ran, but the mageck count step didnt pass the parameter in as far as i can tell?

/usr/local/bin/mageck count -l gRNA_annotated_geneSymbol_withID.tsv -n count_table --sample-label treated_5hr_1,treated_5hr_2,treated_5hr_3,treated_7hr_1,treated_7hr_2,treated_7hr_3,control_1,control_2,control_3 --fastq 20uM-etoposide-5hrs-1_S1_R1_001.fastq.gz 20uM-etoposide-5hrs-2_S4_R1_001.fastq.gz 20uM-etoposide-5hrs-3_S6_R1_001.fastq.gz 20uM-etoposide-7hrs-1_S2_R1_001.fastq.gz 20uM-etoposide-7hrs-2_S2_R1_001.fastq.gz 20uM-etoposide-7hrs-3_S7_R1_001.fastq.gz Untreated-1_S1_R1_001.fastq.gz Untreated-2_S3_R1_001.fastq.gz Untreated-3_S5_R1_001.fastq.gz

oh i was a bit tired and added it for mageck mle 😞 if you’re interested in contributing you could look at my PR and do the same for mageck count? Otherwise i will look at it another day

do you agree it should be on both?

I’m a bit surprised to see it in mageck count actually because for mageck mle it makes the model semi supervised and in my experience then adjusts the pvalues

does it change the normalization on

count

? im not familiar with any of this so will defer to your expertise

--norm-method {none,median,total,control}	Method for normalization, including "none" (no normalization), "median" (median normalization, default), "total" (normalization by total read counts), "control" (normalization by control sgRNAs specified by the --control-sgrna option).

https://sourceforge.net/p/mageck/wiki/usage/#count Looks like it does but i wonder if then the param --norm-method should then be added

i wish their documentation were better 😐

ah that’d be great

otherwise i generally agree that if the negatives are getting passed in, the norm method should also change to match

That's something that we'd be able to set in crisprseq fairly easily

hey @Laurence Kuhlburger so it sounds like there is likely value in allowing the user to turn on the non-targeting gRNA normalization and correction in mle, but it seems probably not good practice to enable those by default

Yes right now I'd add it as a feature