Dear all,
I am currently trying to evaluate UMI recovery, but I have not worked with UMI analysis before.
In my data, some UMIs are extremely dominant, while many others have only 1–2 reads. In this analysis, we are not capturing transcripts, so there is no mapping step. Is this pattern normal, or should I clean or filter the data using a specific protocol?
Do you typically keep UMIs with very low counts, such as 1–5 reads, as well as UMIs with very high counts, for example more than 10,000 reads?
Do you have any suggestion on how to evaluate the UMI in this data?
There can be several reasons why you get such uneven count distribituions.
If have a pool of say X molecules, each at the same concentration, and then draw Y random samples from this pool, with replacement, then the count for each molecule is going to be binomially distributed, which can be approximated by a poisson distribution when the numbers are big enough. Poisson distributions are definately have most things have a small count and a few things with a large count. The extent to which this is the case depends on the ratio of X to Y. If you are seeing most things having only one count, then that probably means that you are under-sequencing, but if you are seeing, say 5 counts on average, then you will probably not gain more by sequencing more.
PCR bias: if you have applied PCR in your library prep, then some moleulces will be amplified more than others. A molecule that gets an early boost in PCR amplifcation is likely to blow up in later cycles as it competes for limiting polymerase and primer. Add to this that the efficiency of amplication in PCR is biased by sequence, particularly at the start of the sequence (where the UMI is) and some sequences will be more likely to be amplified than others. This is infact the primary reason we use UMIs in the first place: the count of reads after PCR is not a very accurate estimate of the number of reads before PCR.
PCR and sequencing errors: PCR and sequencing are not perfect, they can introduce errors. At each cycle of PCR, there is some probability that an error will be introduced into the UMI sequence. If you only have a small number of molecules at the start, but at each cycle of PCR there is a chance of one copy of each UMI mutating, then you could end up with a sitatuion with a small number of high count UMIs and a large number of low count UMIs. We specifically wrote UMI-tools to deal with this.
Of these, really only 3 & 4 are a genuine worry for the quality of the results, and 3 can be solved by using a error aware UMI grouper. While the scripts in the UMI-tools pacakge require reads to be mapped, the package does expose the API so that you can utilise the error correction algorithm directly. see https://umi-tools.readthedocs.io/en/latest/API.html
Thank you for your insight. I initially had 28k UMIs, but after running the clusterer from UMI-tools, only 2.8k UMI clusters remained.
Do you know if this indicate that the true number of unique molecules is much lower than expected, or that the sequencing depth was insufficient and more sequencing is needed?
28K UMIs to 2.8K UMIs is a bigger drop than I would usually expect.
There are two possible explainations.
1) You really do only have 2.8K molecules, and far from having insufficient sequencing depth, you've massively amplified and then over-sequenced.
or
2) You've exhausted the available UMI space, and are collapsing UMIs that shouldn't be collapsed. If you collapse any two UMIs that differ by only a single nucleotide, then you soon start collapsing UMIs that shouldn't be collapsed because they just have similar UMIs by chance. UMI-tools gets around this by only collapsing UMI A into UMI B if the counts for A are much lower than the counts for B. However, even this only withstands overcollapsing for so long. We find that once you've started using more than about 30% of all available UMIs, you start to overcollapse. How long are your UMIs?
For the directional method, I obtained 2.8k UMIs, whereas for the adjacency method I obtained 28.1k UMIs out of a total of 28.5k UMIs. I am reading your paper to understand the difference between these two methods and which one is more appropriate for my case, where UMIs cannot be clustered based on shared genomic coordinates.
Do you think the clusterer function is sufficient for grouping UMIs when there are no shared genomic coordinates?
My UMI is 8bp.
How strange. I would normally expect 'adjacency' to leave you with fewer clusters than 'directional'. I'll have to have a hard think about what sort of situation would leave to it giving you more.
You 8bp UMI gives you 65,536 possible UMIs, and you are seeing 28,500. That suggests to me that you are seeing too many UMIs given the given length. Do you have any other information that might distringuish the reads? Do they all have identical sequences other than the UMI?
oh really, I have no experience with checking UMI before, I thought that 28,500 is quite low compared to expected 65,536. What number of unique UMI that you can expected from 8bp?
below I can quickly show you some code:
cat c_umi.txt | head -10
166666 TGATAATT
100576 AACCGTTG
74428 CTGCATTT
52237 AATCTTTG
50909 CATAGCTA
42313 TCCATTTC
40215 GAAGTCCT
39021 GTCTTTAA
37549 TCGTTCGG
29821 CCATGTTA
cat c_umi.txt | tail -10
1 AAAAAGGT
1 AAAAAGGG
1 AAAAAGGC
1 AAAAAGGA
1 AAAAAGCT
1 AAAAAGAT
1 AAAAAGAG
1 AAAAACGT
1 AAAAACGG
1 AAAAAACG
cat c_umi.txt | head -8000 | tail -10
3 ATTATTAA
3 ATTATGTG
3 ATTATGGA
3 ATTATCTA
3 ATTATCCT
3 ATTATCAA
3 ATTATATT
3 ATTATAGT
3 ATTATAAT
3 ATTATAAC
here with umi_tools:
>>> di_clusterer = UMIClusterer(cluster_method="directional")
>>> adj_clusterer = UMIClusterer(cluster_method="adjacency")
>>>
>>> umis = {}
>>> with open("c_umi.txt") as f:
... for line in f:
... count, seq = line.split()
... umis[seq.encode()] = int(count)
...
...
>>> di_clusterer_umi = di_clusterer(umis, threshold=1)
>>> print(len(di_clusterer_umi))
2132
>>> adj_clusterer_umi = adj_clusterer(umis, threshold=1)
print(len(umis.keys()))
>>> print(len(adj_clusterer_umi))
28145
>>>
>>> print(len(umis.keys()))
28517
What number of unique UMI that you can expected from 8bp?
The number of UMIs you can expect depends on the biology of the system and the chemstriy of the library prep protocol, not the length of the UMI. However, for us to be confident that two reads with the same UMI are PCR duplicates we have to assume that the chance of two biologically independent molecules being attached to the same UMI is low.
If there are only a small number of molecules and a large number of possible UMI sequences, then this is clearly true. Its also clearly true that if you have 100k molecules, and only 65k UMI sequences, some of the UMI sequences have to be used twice.
However, its worse than this becauase if you collapse UMIs that are only 1bp different from each other, then its not just that we have to assume that any two molecules with the same UMI are PCR duplicates, but that two reads with similar UMIs are PCR duplciates.
Our emprical studies suggest that things start to go wrong when you've used around a third of all the possible UMIs. For 65,536 possible UMIs, this means about 21k UMIs is the limit at which I think the assumptions of UMI deduplication start to break down.
Just to be clear, this is usage of UMIs for one particular sequence. If you have different sequences, it would make sense to process each sequence seperately. This is one of the main reasons that UMI-tools considers mapping position as well as UMI, but read sequence could well also work.
Log in to answer this question.
What kind of data is this?
lineage tracing with cas9, we try to capture specific genomic region contain gRNA
There is an extensive answer below that may cover all use cases but it may be interesting to find out where the UMI's are coming from? (adapters for library prep or some other way).