PacBio with 10X barcodes, and partial alignments within it (ideally aiming for all_contig.fasta-like output)

[ktpolanski]

First off, thank you for the package.

There are protocols which mix 10X barcodes with PacBio long read sequencing. The paper uses MiXCR for processing, but more recent versions of the package outperform their proposed workflow in terms of tidiness.

I identify appropriate adapters at the start/end of the read, doing a reverse complement if needed, pull out the cell barcode and UMI from the start of the read into R1.fastq, trim the poly-T tail, reverse complement and export the resulting sequence into R2.fastq. I was able to get the following to clear:

mixcr analyze generic-pacbio-with-umi \
    --species hsa \
    --tag-pattern "^(UMI:N{28})\^(R2:*)" \
    R1.fastq \
    R2.fastq \
    mixcr/mixcr

mixcr exportClones \
    --dont-split-files \
    -isProductive CDR3 \
    -uniqueTagCount Molecule \
    -tagCounts \
    mixcr/mixcr.clns \
    mixcr/mixcr.tsv

This results in a file which in principle holds all required information for importing the data into Scirpy via the AirrCell() route, so there's something. However, the identified contigs are the high quality ones, while sometimes partial alignments are of interest. Ultimately I'd like to get an output akin to all_contig.fasta from Cell Ranger, with reconstructed contigs per cell, including the partial stuff that currently doesn't make it. I've ran into a number of snags while attempting to get this to work:

1. --tag-pattern "^(CELL:N{16})(UMI:N{12})\^(R2:*)" crashes at the refining tags stage. I've found the various yamls and I see no mention of cell in the config, which is understandable as it was calibrated for UMIs. How does one include a cell tag in this?

Please copy the following information along with the stacktrace:
   Version: 4.6.0; built=Sat Dec 09 19:48:42 UTC 2023; rev=c9fafa41fe; lib=repseqio.v4.0
        OS: Linux
      Java: 17.0.6
  Cmd args: refineTagsAndSort --report pls/mixcr.refine.report.txt --json-report pls/mixcr.refine.report.json pls/mixcr.vdjca pls/mixcr.refined.vdjca
picocli.CommandLine$ExecutionException: Error while running command refineTagsAndSort java.lang.IllegalArgumentException: Unexpected filter tag request. Requested UMI but next level is CELL.
        at com.milaboratory.mixcr.cli.Main.registerExceptionHandlers$lambda-12(SourceFile:395)
        at picocli.CommandLine.execute(CommandLine.java:2088)
        at com.milaboratory.mixcr.cli.CommandAnalyze$Cmd$PlanBuilder.executeSteps(SourceFile:463)
        at com.milaboratory.mixcr.cli.CommandAnalyze$Cmd.run0(SourceFile:421)
        at com.milaboratory.mixcr.cli.MiXCRCommand.run(SourceFile:37)
        at picocli.CommandLine.executeUserObject(CommandLine.java:1939)
        at picocli.CommandLine.access$1300(CommandLine.java:145)
        at picocli.CommandLine$RunLast.executeUserObjectOfLastSubcommandWithSameParent(CommandLine.java:2358)
        at picocli.CommandLine$RunLast.handle(CommandLine.java:2352)
        at picocli.CommandLine$RunLast.handle(CommandLine.java:2314)
        at picocli.CommandLine$AbstractParseResultHandler.execute(CommandLine.java:2179)
        at picocli.CommandLine$RunLast.execute(CommandLine.java:2316)
        at com.milaboratory.mixcr.cli.Main.registerLogger$lambda-27(SourceFile:514)
        at picocli.CommandLine.execute(CommandLine.java:2078)
        at com.milaboratory.mixcr.cli.Main.main(SourceFile:101)
Caused by: java.lang.IllegalArgumentException: Unexpected filter tag request. Requested UMI but next level is CELL.
        at com.milaboratory.o.gf$b.getKeyExtractor(SourceFile:957)
        at com.milaboratory.mitool.refinement.gfilter.KeyedFilterContext.getKeyExtractor(SourceFile:97)
        at com.milaboratory.mitool.refinement.gfilter.GroupFilter.filter(SourceFile:297)
        at com.milaboratory.o.gf.c(SourceFile:660)
        at com.milaboratory.mixcr.cli.CommandRefineTagsAndSort$Cmd.run1(SourceFile:353)
        at com.milaboratory.mixcr.cli.MiXCRCommandWithOutputs.run0(SourceFile:69)
        at com.milaboratory.mixcr.cli.MiXCRCommand.run(SourceFile:37)
        at picocli.CommandLine.executeUserObject(CommandLine.java:1939)
        at picocli.CommandLine.access$1300(CommandLine.java:145)
        at picocli.CommandLine$RunLast.executeUserObjectOfLastSubcommandWithSameParent(CommandLine.java:2358)
        at picocli.CommandLine$RunLast.handle(CommandLine.java:2352)
        at picocli.CommandLine$RunLast.handle(CommandLine.java:2314)
        at picocli.CommandLine$AbstractParseResultHandler.execute(CommandLine.java:2179)
        at picocli.CommandLine$RunLast.execute(CommandLine.java:2316)
        at com.milaboratory.mixcr.cli.Main.registerLogger$lambda-27(SourceFile:514)
        at picocli.CommandLine.execute(CommandLine.java:2078)
        ... 13 more

2. While looking at the docs for align, I was able to find the parameters -OallowPartialAlignments=true -OallowNoCDR3PartAlignments=true -OallowChimeras=true. Passing them to the isolated align call has the successfully aligned read count go up by about 50%, but by the time it's time to assemble the actual used read count goes back down to what it was before these parameters were added and the identified clonotype count is similar. How does one have these partial/chimeric reads make it to the output somehow in a collapsed fashion?
3. Is there some way to generate the aforementioned 10X-like contig sequence per cell, ideally with UMI counts to boot?
4. Is there some way to call clonotypes jointly across multiple samples? It seems that MiXCR applies some sophisticated algorithms to the process, and it seems like running samples individually may end up with slightly disjoint clonotype definitions as a result. I tried running multiple samples simultaneously via file name comprehension, but at the end the clonotypes were called independently for each sample. From my understanding, postanalysis overlap doesn't actually link up the clonotypes across samples.

Thanks a lot and sorry for the trouble!

[mizraelson]

Hi,

Can you share the initial library structure? You mentioned identifying the reads by adapters, so I assume you know the sequence that comes before the CELL barcode in the original reads. I can make a dedicated preset that will not require any preprocessing.

Additionally, MiXCR uses partial alignments (those that do not cover the CDR3 region) to build complete clones. The minimal requirement for a clone is that it covers at least the CDR3 region, so the final clonotype table will not include anything that lacks CDR3. It is possible to align such reads and export them as alignments in the early stages of the analysis, but they will not pass clonotype assembly.

Also, I am not sure I understand what you mean by "call clonotypes jointly across multiple samples." Could you elaborate on this?

If you write to us at [email protected], we can arrange a call to discuss your dataset in detail. Let me know if you are interested.

[ktpolanski]

Thank you for your response. Sorry for asking about so much stuff at once.

• A global config sounds quite tempting, but I think the one with a cell barcode present will be sufficient for the time being as I'm processing the reads this way for TRUST4 anyway. Plus the presence or lack thereof of the TSO at the end of the read has become an interesting discussion point with a coworker, which opens various cans of worms. Might be worth doing later but needs some details ironed out and considered on my end, meanwhile the cell one would help immediately.
• Is there some way to convert the alignments to FASTA sequences of the contigs, with cell barcode in tow along with UMI count?
• When analysing single-cell data, it's pretty common to identify clonotypes jointly across multiple samples and examine overlaps, with an example being the repertoire overlap section from the scirpy tutorial. In the interest of internal cohesion, it seems like it would be beneficial to be able to have MiXCR perform its clonotype construction across multiple samples at once. The output can then be ferried into scirpy via the AirrCell() route for downstream analysis.

[mizraelson]

Hi,

1. What output do you expect to get with TRUST4 compared to using MiXCR? I would strongly recommend making the preset for the intact data. In any case to analyze this data you will need a custom preset (attached). From the tag pattern I assume you used the 10x V3 kit? To use the preset put the YAML file into the directory you run MiXCR from or to ~/.mixcr/presets/ folder.
   Then run:

mixcr analyze local:10x-pac-bio-split-v1 \
    --species hsa \
    input_R1.fastq.gz \
    input_R2.fastq.gz \
    output_prefix

2. The alignments can't be exported as a FASTA file. You can export them in a tabular format using:

mixcr exportAlignments \
    input.vdjca \
    alignments.tsv

In this file you will have an alignment for every read with a CELL and UMI barcode sequence. The alignments are not aggregated clones yet, so you can't see the UMI count here.

The custom preset by default will produce the clonotype table with a CELL barcode and a UMI count for every assembled clone.

3. MiXCR has a couple of ways to make an overlap. There is no need to process multiple samples at once. You can use clns files from multiple samples to overlap the repertoires. The mixcr exportClonesOverlap command overlaps multiple samples by a certain criteria (e.g. CDR3 sequence) and returns a table of occurrences for each sample. Also you can use mixcr postanalysis overlap to compare the repertoires using different metrics like Pearson's correlation etc.

preset.zip

[ktpolanski]

Can confirm the preset worked like a charm, thank you very much! 4.6.0 did not recognise assembleCells, but I downloaded the development version and it all worked fine.

In the docs, the assemble step stresses that it performs a lot of error correction. As such, I'm concerned that running this separately on multiple samples may lead to different local optima in some cases, and processing the data jointly would have a single clonotype inferred. I guess I can fix it in post by doing a thresholded Hamming distance if things refuse to cooperate.

As for TRUST4, it creates a fasta representation of the contigs, including ones that only have the J+C genes. Which is the output I desire.

[mizraelson]

Yes, good thinking with using the development version!

The assemble step indeed performs several levels of correction: UMI-based error correction, sequencing error correction (mapping), and PCR error correction (clustering). However, you can't cluster different samples together; otherwise, you will end up mixing the clones. It is possible to pass multiple FASTQ files to MiXCR and process them as a single sample, but then you won't be able to tell which clone came from which file because they would be mixed, and a single clone would be assembled from reads from multiple files. In other words, if you want to compare the repertoires, they should be processed independently and then compared. If you want to cluster clones between different repertoires, that is also possible, but there would not be different samples anymore.

We have a dedicated software, MiTool, which you can use to export consensuses as a FASTQ file (but not aligned to the reference of course).

Additionally, if you use the --keep-non-CDR3-alignments parameter with the analyze command, the alignments will contain partially covered sequences too (C+J). You will see them in the output of the exportAlignments command.

[ktpolanski]

I gave exportAlignments a try, using the following command:

mixcr exportAlignments -uniqueTagCount Molecule mixcr.refined.vdjca mixcr.refined.tsv

This created a table which includes a prominent sequence in the targetSequences column, which I can easily peel out without further software, but seems to come with a singleton tagValueCELL and tagValueUMI per row, further validated by uniqueMoleculeCount having a value of 1 all the way through. This is consistent with what you wrote earlier, that this summary is on a read level. So this would still need some manner of cross-UMI collapsing on a cell level to get the desired output, which you do in assemble for the complete ones.

It's fine if this is out of scope! The custom preset helps plenty already.

Is MiTool this thing? Seeing how its last release is from 2017, I had to conda create --name oldjdk -c conda-forge openjdk=7.0.161 to get the output of mitools -h without it complaining about the java version. Once successful, I see nothing there that seems to relate to what we're discussing.

[mizraelson]

Yes, we have an updated version of MiTool, but it is not public. Please write to [email protected] and I will share the instructions.

[ktpolanski]

Sorry to bump this, but we've gotten some fresh data and I'm running the preset you so nicely shared. I left things chugging, came back a day later expecting to find the entire lot of samples completed, but instead found the second sample stuck printing "Initializing: 100%" over and over again in the barcode correction step rather than moving on to anything. I tried skipping the second sample, and then the third one got stuck in a similar fashion, though I didn't wait a full day this time. We've got more reads this time around, could this somehow be related?

Ran like so for transparency:

/mnt/mixcr-develop/mixcr analyze local:10x-pac-bio-split-v1 \
    --species hsa \
    --set-whitelist CELL=file:~/STAR-2.7.3a/whitelists/arc.txt \
    R1.fastq \
    R2.fastq \
    mixcr-bespoke-preset/mixcr

[ktpolanski]

I'm pretty sure I won't be able to get clearance, sorry. I'm running the same MiXCR as I obtained last time, the July development version, in case it matters.

[mizraelson]

Can you try adding both: --dont-correct-tag-type Cell and --dont-correct-tag-type Molecule. And lets wait for a while. Also, just in case, you can try updating to the latest develop version.

[ktpolanski]

Sure, I'll poke around with this and return with findings. Strange stuff going on!

Also don't fret about this too much, worst comes to worst the thing that falls out of my initial command can be post-processed into useable output, so we're not completely stuck.

[mizraelson]

Yes, but if you use ^(UMI:N{28})\^(R2:*) you treat each UMI as a unique cell. So you will end up with way more cells than there is instead of making consensuses within each cell.

[ktpolanski]

import pandas as pd
import numpy as np

#get a cell by clonotype count matrix

df = pd.read_table("mixcr-pacbio-umi/mixcr.tsv", index_col=0)
df.index = [str(i) for i in df.index]

ccm = {}
for ind in df.index:
        #the records are formatted in one of these two ways:
        #{ACTG=2.0,ACGT=1.0}
        #ATGC: 1.0
        #can get one entry per "UMI" (CB+UB) by splitting on ,
        #then can strip out all sorts of junk by splitting on =, : and lstripping {
        #and that's it, that's the "UMI". take the first 16bp for CB
        #then need to pd.Series().value_counts() it to get counts per CB
        #and then turn it to a dict for ease of df creation later
        ccm[ind] = pd.Series([i.split("=")[0].split(":")[0].lstrip("{")[:16] for i in df.loc[ind, 'tagCounts'].split(",")]).value_counts().to_dict()

ccm_df = pd.DataFrame.from_dict(ccm).fillna(0)