generic-amplicon-with-umi has a lot less clone than the rna-seq pipeline

[li1311139481]

Recently I have a batch of experimental data, 5' contains UMI. I used two preset processes, rna-seq and generic-amplicon-with umi, to analyze the data. The comparison results are as follows

First question:

Firstly, I use rna-seq pipeline. Because of the existence of umi structure, umi was first removed by cutadapt, and then analyzed by rna-seq pipeline

mixcr analyze rna-seq --species mmu --threads 8 --tag-pattern "^(UMI:N{10})tttcttatatggg\^(R2:*)" --use-local-temp ${BASE}_1_trm.fq.gz ${BASE}_2_trm.fq.gz /work/ZF091_092_093/mixcr/$BASENAME

another, I use generic-amplicon-with-umi pipeline.
mixcr analyze generic-amplicon-with-umi --species mmu --threads 8 --rna --rigid-left-alignment-boundary --floating-right-alignment-boundary C --use-local-temp --tag-pattern '^(UMI:N{10})tttcttatatggg(R1:*)\^(R2:*)' ${BASE}_1.fq.gz ${BASE}_2.fq.gz /work/ZF091_092_093/mixcr_deumi/$BASENAME

Then i will show the report

1. For rna-seq pipeline

ZF091_CD62L_a_i502_split.align.report.txt	Successfully aligned reads: 4335007 (92.2%)
ZF091_CD62L_b_i502_split.align.report.txt	Successfully aligned reads: 4406064 (92.75%)
ZF091_KLRG1_a_i502_split.align.report.txt	Successfully aligned reads: 4891512 (92.62%)
ZF091_KLRG1_b_i502_split.align.report.txt	Successfully aligned reads: 5482410 (93.27%)
ZF092_pool_a_Colon1_split.align.report.txt	Successfully aligned reads: 1448789 (93.19%)
ZF092_pool_a_Colon2_split.align.report.txt	Successfully aligned reads: 1281754 (93.54%)
ZF092_pool_a_Tumor1_split.align.report.txt	Successfully aligned reads: 1194369 (92.66%)
ZF092_pool_a_Tumor2_split.align.report.txt	Successfully aligned reads: 84721 (93.74%)
ZF092_pool_b_Colon1_split.align.report.txt	Successfully aligned reads: 1247662 (92.95%)
ZF092_pool_b_Colon2_split.align.report.txt	Successfully aligned reads: 1364218 (92.91%)
ZF092_pool_b_Tumor1_split.align.report.txt	Successfully aligned reads: 1296362 (93.75%)
ZF092_pool_b_Tumor2_split.align.report.txt	Successfully aligned reads: 193452 (93.31%)
ZF093_pool_a_a_Ap_split.align.report.txt	Successfully aligned reads: 1724058 (94.31%)
ZF093_pool_a_a_Apt_split.align.report.txt	Successfully aligned reads: 754825 (95.6%)
ZF093_pool_a_a_Tp_split.align.report.txt	Successfully aligned reads: 694385 (94.14%)
ZF093_pool_a_a_Tpt_split.align.report.txt	Successfully aligned reads: 1290968 (95.97%)
ZF093_pool_b_b_Ap_split.align.report.txt	Successfully aligned reads: 1977512 (92.8%)
ZF093_pool_b_b_Apt_split.align.report.txt	Successfully aligned reads: 678650 (92.94%)
ZF093_pool_b_b_Tp_split.align.report.txt	Successfully aligned reads: 652086 (93.1%)
ZF093_pool_b_b_Tpt_split.align.report.txt	Successfully aligned reads: 742757 (92.89%)

2. For generic-amplicon-with-umi pipeline

ZF091_CD62L_a_i502_split.align.report.txt	Successfully aligned reads: 3898426 (82.91%)
ZF091_CD62L_b_i502_split.align.report.txt	Successfully aligned reads: 3333212 (70.16%)
ZF091_KLRG1_a_i502_split.align.report.txt	Successfully aligned reads: 4343059 (82.22%)
ZF091_KLRG1_b_i502_split.align.report.txt	Successfully aligned reads: 3905528 (66.43%)
ZF092_pool_a_Colon1_split.align.report.txt	Successfully aligned reads: 1295971 (83.36%)
ZF092_pool_a_Colon2_split.align.report.txt	Successfully aligned reads: 1177485 (85.93%)
ZF092_pool_a_Tumor1_split.align.report.txt	Successfully aligned reads: 1128498 (87.55%)
ZF092_pool_a_Tumor2_split.align.report.txt	Successfully aligned reads: 70126 (77.59%)
ZF092_pool_b_Colon1_split.align.report.txt	Successfully aligned reads: 890395 (66.33%)
ZF092_pool_b_Colon2_split.align.report.txt	Successfully aligned reads: 1021108 (69.53%)
ZF092_pool_b_Tumor1_split.align.report.txt	Successfully aligned reads: 1071491 (77.48%)
ZF092_pool_b_Tumor2_split.align.report.txt	Successfully aligned reads: 174505 (84.17%)
ZF093_pool_a_a_Ap_split.align.report.txt	Successfully aligned reads: 1511325 (82.67%)
ZF093_pool_a_a_Apt_split.align.report.txt	Successfully aligned reads: 658170 (83.35%)
ZF093_pool_a_a_Tp_split.align.report.txt	Successfully aligned reads: 607537 (82.37%)
ZF093_pool_a_a_Tpt_split.align.report.txt	Successfully aligned reads: 1183112 (87.94%)
ZF093_pool_b_b_Ap_split.align.report.txt	Successfully aligned reads: 1430138 (67.11%)
ZF093_pool_b_b_Apt_split.align.report.txt	Successfully aligned reads: 510113 (69.85%)
ZF093_pool_b_b_Tp_split.align.report.txt	Successfully aligned reads: 472245 (67.42%)
ZF093_pool_b_b_Tpt_split.align.report.txt	Successfully aligned reads: 531577 (66.47%)

Result
Although input data for the rna-seq pipeline were pre-processed by cutadapt, the proportion of meaningful reads was higher. However, the written ratio of my cutadapt is as follows: the mapping ratio is 0.91*0.92 ~=83% according to the ratio. However, the mapping ratio of generic-amplicon-with-umi pipeline has been as low as 66.33%

ZF091_CD62L_a_i502_split.cutadapt.log	Total written (filtered):  1,185,183,363 bp (91.5%)
ZF091_CD62L_b_i502_split.cutadapt.log	Total written (filtered):  1,196,348,032 bp (91.4%)
ZF091_KLRG1_a_i502_split.cutadapt.log	Total written (filtered):  1,328,010,368 bp (91.4%)
ZF091_KLRG1_b_i502_split.cutadapt.log	Total written (filtered):  1,478,614,260 bp (91.3%)
ZF092_pool_a_Colon1_split.cutadapt.log	Total written (filtered):    392,378,402 bp (91.5%)
ZF092_pool_a_Colon2_split.cutadapt.log	Total written (filtered):    345,688,013 bp (91.5%)
ZF092_pool_a_Tumor1_split.cutadapt.log	Total written (filtered):    325,237,033 bp (91.5%)
ZF092_pool_a_Tumor2_split.cutadapt.log	Total written (filtered):     22,819,147 bp (91.6%)
ZF092_pool_b_Colon1_split.cutadapt.log	Total written (filtered):    338,374,019 bp (91.5%)
ZF092_pool_b_Colon2_split.cutadapt.log	Total written (filtered):    369,951,706 bp (91.5%)
ZF092_pool_b_Tumor1_split.cutadapt.log	Total written (filtered):    348,472,735 bp (91.4%)
ZF092_pool_b_Tumor2_split.cutadapt.log	Total written (filtered):     52,296,503 bp (91.5%)
ZF093_pool_a_a_Ap_split.cutadapt.log	Total written (filtered):    461,485,374 bp (91.6%)
ZF093_pool_a_a_Apt_split.cutadapt.log	Total written (filtered):    199,493,300 bp (91.6%)
ZF093_pool_a_a_Tp_split.cutadapt.log	Total written (filtered):    186,220,631 bp (91.5%)
ZF093_pool_a_a_Tpt_split.cutadapt.log	Total written (filtered):    340,019,327 bp (91.6%)
ZF093_pool_b_b_Ap_split.cutadapt.log	Total written (filtered):    537,630,686 bp (91.5%)
ZF093_pool_b_b_Apt_split.cutadapt.log	Total written (filtered):    184,061,899 bp (91.5%)
ZF093_pool_b_b_Tp_split.cutadapt.log	Total written (filtered):    176,402,082 bp (91.4%)
ZF093_pool_b_b_Tpt_split.cutadapt.log	Total written (filtered):    201,379,798 bp (91.5%)

Second question:

After summarizing the results of the two pipelines, I found some problems.
Summary table is as follows:

cell_number: the number of cells in the sample
lines number of result from rna-seq: The number of rows in the output table of mixcr from rna-seq pipeline
lines number of result from generic-amplicon-with-umi : The number of rows in the output table of mixcr from generic-amplicon-with-umi
unique_CDR3_number from rna-seq: The amino acid sequences of the aaSeqCDR3 columns were de-duplicated, and then how many unique sequences were counted. Results from the rna-seq pipeline
unique_CDR3_number from generic-amplicon-with-umi: The amino acid sequences of the aaSeqCDR3 columns were de-duplicated, and then how many unique sequences were counted. Results from the generic-amplicon-with-umi pipeline.
sum of uniqueMoleculeCount from generic-amplicon-with-umi: Sum the uniqueMoleculeCount column to calculate how many UMIs were detected in total

<style> </style>

name	cell_number	lines number of result from rna-seq	lines number of result from generic-amplicon-with-umi	unique_CDR3_number from rna-seq	unique_CDR3_number from generic-amplicon-with-umi	sum of uniqueMoleculeCount from generic-amplicon-with-umi
ZF091_CD62L_a_i502_deumi_split.clones_TRA.tsv	10000	8087	2303	7424	2219	2682
ZF091_CD62L_b_i502_deumi_split.clones_TRB.tsv	10000	12019	3946	11054	3842	4795
ZF091_KLRG1_a_i502_deumi_split.clones_TRA.tsv	10000	5801	1559	5314	1448	4874
ZF091_KLRG1_b_i502_deumi_split.clones_TRB.tsv	10000	6353	1903	5789	1781	7149
ZF092_pool_a_Colon1_deumi_split.clones_TRA.tsv	10000	4627	938	4264	911	1235
ZF092_pool_a_Colon2_deumi_split.clones_TRA.tsv	10000	4973	1450	4591	1398	1641
ZF092_pool_a_Tumor1_deumi_split.clones_TRA.tsv	10000	3927	305	3634	302	425
ZF092_pool_a_Tumor2_deumi_split.clones_TRA.tsv	5000	1913	383	1812	378	550
ZF092_pool_b_Colon1_deumi_split.clones_TRB.tsv	10000	5527	746	5218	736	1033
ZF092_pool_b_Colon2_deumi_split.clones_TRB.tsv	10000	7351	2472	6884	2438	2885
ZF092_pool_b_Tumor1_deumi_split.clones_TRB.tsv	10000	5082	307	4803	305	419
ZF092_pool_b_Tumor2_deumi_split.clones_TRB.tsv	5000	2662	349	2570	347	469
ZF093_pool_a_a_Ap_deumi_split.clones_TRA.tsv	10000	1355	390	1256	379	778
ZF093_pool_a_a_Apt_deumi_split.clones_TRA.tsv	2000	528	77	495	75	231
ZF093_pool_a_a_Tp_deumi_split.clones_TRA.tsv	900	1274	42	1221	42	75
ZF093_pool_a_a_Tpt_deumi_split.clones_TRA.tsv	500	710	152	678	151	549
ZF093_pool_b_b_Ap_deumi_split.clones_TRB.tsv	10000	2439	654	2257	637	1745
ZF093_pool_b_b_Apt_deumi_split.clones_TRB.tsv	2000	817	92	769	87	422
ZF093_pool_b_b_Tp_deumi_split.clones_TRB.tsv	900	1467	54	1415	54	92
ZF093_pool_b_b_Tpt_deumi_split.clones_TRB.tsv	500	1021	133	978	130	447

1. Compared with rna-seq pipeline, generic-amplicon-with-umi seems to get fewer result rows. Does generic-amplicon-with-umi lose more potential CDR3? By comparing the number of unique_CDR3_number, I found that generic-amplicon-with-umi obtained less than one-third of the amount of aaSeqCDR3 as rna-seq.
   Through manual comparison, I found that in the rna-seq pipeline, as the number of cloneCount in ran-seq decreased, the proportion of corresponding aaSeqCDR3 that could not be searched in generic-amplicon-with-umi increased. If these clonecounts are all from the same mRNA and have the same UMI, I should at least see uniqueMoleculeCount=1 in the generic-amplicon-with-umi result.
   From the above, can I conclude that the generic-amplicon-with-umi process loses more CDR3, especially for clone with a low count? Can you explain the reasons for this? Or if you can help me modify the code to increase the proportion of CDR3 detected and not discard clones with low counts.

2. If an UMI represents a TCR mRNA. Based on the "sum of uniqueMoleculeCount from generic-amplicon-with-umi" column in the table, do you think I am capturing TCR mRNA too inefficiently compared to the number of cells in the sample? Is this because of my experimental procedures?

[mizraelson]

The crucial aspect to consider is the type of data you possess. Is it fragmented (shotgun) bulk RNA-seq, or is it amplicon data (e.g., 5'RACE) with UMI? The presets for these two protocols differ significantly in several ways. To summarize:

• The aligners used are different.
• The rna-seq preset does not process UMIs, even if you specify them within the pattern.
• The rna-seq preset includes additional steps that are only relevant for fragmented data.
• The generic-amplicon-with-umi preset not only corrects data using UMIs but also filters out artificial diversity based on the reads/UMI density plot.
• The rna-seq preset initially assembles clones by CDR3 and then attempts to extend the alignment to cover as much of the receptor sequence as possible. In contrast, the generic-amplicon-with-umi preset, by default, assembles clones only by CDR3.

It's important to note that one of the primary challenges in TCR/BCR repertoire analysis is eliminating the artificial diversity that is invariably present in the sample due to PCR and sequencing errors. A higher number of clones does not necessarily indicate that a certain preset is more effective if it was chosen incorrectly with respect to the data structure.

To address your question:

1. generic-amplicon-with-umi better handles artificial diversity if you data is amplicon based and not fragmented. The algorithms preserve even low abundant clones if those are true clones.
2. It is hard to say anything about the efficiency without knowing the details of the experiment. I see you have different cell populations and tissues which might explain the difference in the number of recovered UMIs per sample.

[li1311139481]

Thanks for the summary, I learned a lot.

1. My data should be what you call amplicon data (5' RACE). Because I do contain UMI and template switching oligos at the 5' end. And there was no fragmentation step before final library construction. So I'm going to try the generic-amplicon-with-umi pipeline

2. you said

The aligners used are different

I understand: Because of the different aligners used, there are different comparison rates. Right? I'm really worried that the mapping ratio. I noticed ALERT in the output report, does this mean that the data is not of high quality and I should not put it in the paper

Successfully aligned reads: 3898426 (82.91%)
  Successfully aligned reads:                           82.91% [WARN]
Successfully aligned reads: 3333212 (70.16%)
  Successfully aligned reads:                           70.16% [WARN]
Successfully aligned reads: 4343059 (82.22%)
  Successfully aligned reads:                           82.22% [WARN]
Successfully aligned reads: 3905528 (66.43%)
  Successfully aligned reads:                           66.43% [ALERT]
Successfully aligned reads: 1295971 (83.36%)
  Successfully aligned reads:                           83.36% [WARN]
Successfully aligned reads: 1177485 (85.93%)
  Successfully aligned reads:                           85.93% [OK]
Successfully aligned reads: 1128498 (87.55%)
  Successfully aligned reads:                           87.55% [OK]
Successfully aligned reads: 70126 (77.59%)
  Successfully aligned reads:                           77.59% [WARN]
Successfully aligned reads: 890395 (66.33%)
  Successfully aligned reads:                           66.33% [ALERT]
Successfully aligned reads: 1021108 (69.53%)
  Successfully aligned reads:                           69.53% [ALERT]
Successfully aligned reads: 1071491 (77.48%)
  Successfully aligned reads:                           77.48% [WARN]
Successfully aligned reads: 174505 (84.17%)
  Successfully aligned reads:                           84.17% [WARN]
Successfully aligned reads: 1511325 (82.67%)
  Successfully aligned reads:                           82.67% [WARN]
Successfully aligned reads: 658170 (83.35%)
  Successfully aligned reads:                           83.35% [WARN]
Successfully aligned reads: 607537 (82.37%)
  Successfully aligned reads:                           82.37% [WARN]
Successfully aligned reads: 1183112 (87.94%)
  Successfully aligned reads:                           87.94% [OK]
Successfully aligned reads: 1430138 (67.11%)
  Successfully aligned reads:                           67.11% [ALERT]
Successfully aligned reads: 510113 (69.85%)
  Successfully aligned reads:                           69.85% [ALERT]
Successfully aligned reads: 472245 (67.42%)
  Successfully aligned reads:                           67.42% [ALERT]
Successfully aligned reads: 531577 (66.47%)
  Successfully aligned reads:                           66.47% [ALERT]

2. your said

The rna-seq preset initially assembles clones by CDR3 and then attempts to extend the alignment to cover as much of the receptor sequence as possible. In contrast, the generic-amplicon-with-umi preset, by default, assembles clones only by CDR3.
generic-amplicon-with-umi better handles artificial diversity if you data is amplicon based and not fragmented. The algorithms preserve even low abundant clones if those are true clones.

I understand: Although rna-seq pipeline gets a lot of clones, that is due to the fact that rna-seq pipeline's preset process cannot detect base mutations caused by the PCR and sequencing process, e.g., from the same mRNA molecule, and if a base mutation occurs during the PCR process, the rna-seq pipeline process assembles the mutated sequence into a new clone, but generic-amplicon-with-umi avoids this by using UMI.
If I understand correctly, I think this explains why the results obtained in the rna-seq pipeline, as the cloneCount decreases, cannot be searched for in the output of generic-amplicon-with-umi

Another Question

I found that although my data is 5' RACE data, the cDNA is not synthesized to the 5' end of the TCR V region during reverse transcription of the mRNA. So is it possible for me to change the --rigid-left-alignment-boundary option to --floating-left-alignment-boundary? Will this be able to increase my mapping ratio? Also can you by the way explain what <anchor_point> is under the --floating-left-alignment-boundary parameter? how to specify it ? I really don't understand it. Thank you!

[mizraelson]

1. The alignment rate for most of your samples exceeds 80%, which is great. The samples at 66% is somewhat on the lower side. but it's still perfectly fine to use for your paper (believe me, I have seen way worse data published) because MiXCR effectively recovers reliable information across a wide spectrum of data quality. The ALERT is simply a prompt suggesting you might want to investigate further. In the alignment report, you can find details on why some reads were discarded.

2. The rna-seq preset also uses specialized algorithms to deal with sequencing and PCR errors, but UMIs significantly enhance this efficiency. And, you're right, that using the rna-seq preset tends to retain more noise that you see as more "clones" in the output.

3. Based on your tag pattern '^(UMI:N{10})tttcttatatggg(R1:*)\^(R2:*)', it appears you're trimming both the UMI and the adapter sequences. If no artificial sequences remain, you can use --rigid-left-alignment-boundary. This command instructs MiXCR to align precisely up to the last nucleotide (before tttcttatatggg sequence) on the 5' end. While the <anchor_point> can be adjusted to target different sections of the V gene, the absence of adapter sequences means you can safely stick with the default setting.

[li1311139481]

Thank you most sincerely. I can confidently just use the default parameters of generic-amplicon-with-umi.