Error: max() arg is an empty sequence with amplicon (high coverage panels) data

[naumenko-sa]

Dear CRISPresso2 team!

Thanks for developing the tool and helping me so quickly with the WGS case!

Describe the bug
This time I am running CRISPRessoPooled for high coverage panel (~10000X, ~115 baits 121bp wide) data in Amplicons or Genome mode and
getting Error: max() arg is an empty sequence error.
I am using a docker image and it worked well with WGS data from another experiment.

Expected behavior
Aligning reads to the genome, discovering covered regions, calling editing events.

To reproduce

singularity run \
-e crispresso2_latest.sif \
CRISPRessoPooled \
-r1 ${1}_1.fq -r2 ${1}_2.fq \
-x /path/to/reference/mm10_plus/bowtie2/mm10_plus \
-p 10 -n ${1}_only_genome \
--debug

Debug output

Aligning reads to the provided genome index...
aligning with command: bowtie2 -x /path/to/reference/mm10_plus/bowtie2/mm10_plus -p 10  --end-to-end -N 0 --np 0 --mp 3,2 --score-min L,-5,-1.2  -U CRISPRessoPooled_on_sample_only_genome/out.extendedFrags.fastq.gz 2>>CRISPRessoPooled_on_sample_only_genome/CRISPRessoPooled_RUNNING_LOG.txt| samtools view -bS - | samtools sort -@ 10 - -o CRISPRessoPooled_on_sample_only_genome/sample_only_genome_GENOME_ALIGNED.bam
1681504 reads; of these:
  1681504 (100.00%) were unpaired; of these:
    673 (0.04%) aligned 0 times
    1607636 (95.61%) aligned exactly 1 time
    73195 (4.35%) aligned >1 times
99.96% overall alignment rate
Deleting partially-completed demultiplexing in CRISPRessoPooled_on_sample_only_genome/MAPPED_REGIONS/...
Preparing to demultiplex reads aligned to the genome...
Demultiplexing reads by location (67 genomic regions)...
Parsing the demultiplexed files and extracting locations and reference sequences...
Running CRISPResso on the regions discovered...
Running CRISPResso with 10 processes
Finished all regionsERROR: max() arg is an empty sequence

REPORT_READS_ALIGNED_TO_GENOME_ONLY.txt is empty. In the bam aligned with bwa mem I can see the reads aligned over the targets with >1000X coverage.

I also tried to extract the target sequences (baits) from the reference, split them into small chunks (40bp) and run the Pooled script - same result - none of the regions has enough reads. Reads are 138-140bp in lengths.

Could you please suggest what might be wrong in the setting up the parameters for this run?
Thanks!
Sergey

[kclem]

I believe this error could be coming from one of the sub CRISPResso commands run on each region. You can search for the file that contains this error by changing to the CRISPRessoPooled output directory and using the command grep "max() arg is an empty sequence" CRISPResso_on_*/*CRISPResso_RUNNING_LOG.txt. Does that command produce any output or point you to any files? Ideally I'm looking for a line number to trace this error back to so I can fix it.

[naumenko-sa]

hello @kclem !

Thanks for the quick response!

I'm trying to process these data in 3 modes.

1. Genome mode (-x)
   here all samples fail, there is no sub-folders in CRISPResso_on_sampleX. The only log file is
   CRISPResso_on_sampleX/CRISPRessoPooled_RUNNING_LOG.txt
   here is the tail of it:

Aligning reads to the provided genome index...
aligning with command: bowtie2 -x /path/to/bowtie2/mm10_plus -p 10  --end-to-end -N 0 --np 0 --mp 3,2 --score-min L,-5,-1.2  -U CRISPRessoPooled_on_sample_only_genome/out.extendedFrags.fastq.gz 2>>CRISPRessoPooled_on_sample_only_genome/CRISPRessoPooled_RUNNING_LOG.txt| samtools view -bS - | samtools sort -@ 10 - -o CRISPRessoPooled_on_sample_only_genome/sample_only_genome_GENOME_ALIGNED.bam
1681504 reads; of these:
  1681504 (100.00%) were unpaired; of these:
    673 (0.04%) aligned 0 times
    1607636 (95.61%) aligned exactly 1 time
    73195 (4.35%) aligned >1 times
99.96% overall alignment rate
Deleting partially-completed demultiplexing in CRISPRessoPooled_on_sample_only_genome/MAPPED_REGIONS/...
Preparing to demultiplex reads aligned to the genome...
Demultiplexing reads by location (67 genomic regions)...
Parsing the demultiplexed files and extracting locations and reference sequences...
Running CRISPResso on the regions discovered...
Running CRISPResso with 10 processes
Finished all regionsERROR: max() arg is an empty sequence

There is no dump or stack trace in the log, but there is one in the job output when running with --debug:

INFO  @ Thu, 10 Jun 2021 00:55:20:
         Parsing the demultiplexed files and extracting locations and reference sequences...
INFO  @ Thu, 10 Jun 2021 00:55:20:
         Running CRISPResso on the regions discovered...
INFO  @ Thu, 10 Jun 2021 00:55:20:
         Running CRISPResso with 10 processes
INFO  @ Thu, 10 Jun 2021 00:55:20:
         Finished all regions
CRITICAL @ Thu, 10 Jun 2021 00:55:21:
ERROR: max() arg is an empty sequenceTraceback (most recent call last):
  File "/opt/conda/lib/python2.7/site-packages/CRISPResso2-2.1.1-py2.7-linux-x86_64.egg/CRISPResso2/CRISPRessoPooledCORE.py", line 1316, in main
    CRISPRessoPlot.plot_reads_total(plot_root,df_summary_quantification,save_png,args.min_reads_to_use_region)
  File "/opt/conda/lib/python2.7/site-packages/CRISPResso2-2.1.1-py2.7-linux-x86_64.egg/CRISPResso2/CRISPRessoPlot.py", line 1734, in plot_reads_total
    if max(df['Reads_total'] > 100000):
ValueError: max() arg is an empty sequence

2. Amplicon mode with narrow baits (40bp)
   same as 1, no subfolders,
   error message: 'Skipping the folder : not enough reads, incomplete, or empty folder'

3. Amplicon mode with wide baits (>121bp, baits were sorted and merged if overlapping)

Most samples actually work!

But some samples give another error:

Align reads to the amplicons...
Alignment command: bowtie2 -x CRISPRessoPooled_on_sample_merged_baits/CUSTOM_BOWTIE2_INDEX -p 10  --end-to-end -N 0 --np 0 --mp 3,2 --score-min L,-5,-1.2  -U CRISPRessoPooled_on_sample_merged_baits/out.extendedFrags.fastq.gz 2>>CRISPRessoPooled_on_sample_merged_baits/CRISPRessoPooled_RUNNING_LOG.txt | samtools view -bS - > CRISPRessoPooled_on_sample_merged_baits/CRISPResso_AMPLICONS_ALIGNED.bam
2145972 reads; of these:
  2145972 (100.00%) were unpaired; of these:
    25746 (1.20%) aligned 0 times
    2119038 (98.74%) aligned exactly 1 time
    1188 (0.06%) aligned >1 times
98.80% overall alignment rate
ERROR: list index out of range
samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1

Sergey

[naumenko-sa]

update:
for the samples and probes where CRISPResso2 generated a report, it reports 100% reads modified!
That is surprising, since in the WGS project there is always a fraction < 100%.

[kclem]

Ok. I have fixed the ERROR: max() arg is an empty sequence bug -- but this isn't your main problem. Your main problem is that insufficient reads are aligning to the locations where you think they should be aligning.

CRISPRessoPooled was designed to analyze pooled amplicon sequencing experiments -- is that what your data is? PCR primers for each amplicon create reads that start and stop at exactly the same genomic position. If your reads don't start at the given start position, they will be discarded and you'll end up with an error or empty data as above.

In the genome-only mode (-x) check the produced bam file "CRISPRessoPooled_on_sample_only_genome/sample_only_genome_GENOME_ALIGNED.bam" to make sure that they are aligning to the proper sequences. Then, for one region which has sufficient coverage (>1000x) pull out the genomic region corresponding to one of the reads, and use that for CRISPRessoPooled in amplicon mode (running with this single amplicon). Is this single amplicon sequence one of the ones you provided in the amplicons mode? If not, there may be a problem with the way you define your amplicons. Does this work for you?

[naumenko-sa]

Hi @kclem !

Thanks for the fix and the explanation!
You are right - this not amplicon sequencing, it is rather high coverage panels, i.e. reads don't start at the same position.

the file you mentioned CRISPRessoPooled_on_sample_only_genome/sample_only_genome_GENOME_ALIGNED.bam does have reads aligned to the probes.

So the solution would be just to trim all the reads on the same borders - that would make this data similar to amplicons!

Sergey

[kclem]

Yeah -- you should try CRISPRessoWGS -- this will automatically trim reads that completely overlap the query regions and analyze them using CRISPResso. Just make sure that your query regions are shorter than the read length (otherwise no reads will completely cover the query regions). Good luck!