brwnj / umitools Goto Github PK

View Code? Open in Web Editor NEW

30.0 6.0 8.0 56 KB

Tools to handle reads sequenced with unique molecular identifiers (UMIs).

Home Page: http://brwnj.github.io/umitools

License: MIT License

Python 100.00%

umitools's Introduction

UMI Tools

Tools to handle reads sequenced with unique molecular identifiers (UMIs).

Trim the UMI

Incorporate the UMI into the read name in order to later identify while processing mapped reads.

umitools trim --end 5 unprocessed_fastq NNNNNV > out.fq

If you want to save reads with invalid UMI sequences, you can specify --invalid.

umitools trim --end 5 --invalid bad_umi.fq unprocessed_fastq NNNNNV > out.fq

Remove Duplicates

For any given start site, save only one read per UMI. Writes bed3+ to stdout with before and after counts per start.

umitools rmdup unprocessed.bam out.bam > before_after.bed

Specifying --mismatches will, for a given start site, merge all UMIs within that edit distance into a single unique hit. For example, if a new UMI is within a single mismatch of any existing observed UMIs for a start position, it will be merged and considered a duplicate. The mismatch can occur at any position, regardless of the IUPAC sequence you're using.

Installation

umitools has two requirements: pysam and editdist. Use pip to install pysam.

pip install pysam

editdist has to be downloaded and installed from source (Downloads page).

wget https://py-editdist.googlecode.com/files/py-editdist-0.3.tar.gz
tar xzf py-editdist-0.3.tar.gz
cd py-editdist-0.3/
python setup.py install

Finally download and install umitools from source.

wget -O umitools-master.zip https://github.com/brwnj/umitools/archive/master.zip
unzip umitools-master.zip
cd umitools-master
python setup.py install

umitools's People

Contributors

Stargazers

Watchers

Forkers

jdblischak pombredanne dgennert sunbymoon leangreen scchess geoo95 shulp2211

umitools's Issues

Feature request: Return reads that do not contain a UMI

I would like to investigate the reads that do not contain a UMI. I am curious if they are able to be mapped to the genome. Would it be possible to add a flag to trim.sh to send reads without a UMI to standard error? It would be mutually exclusive with the current --verbose flag.

Here's an idea of how it would be implemented on the command line.

# Print UMI stats to stderr
umitools trim --verbose unprocessed_fastq NNNNNV 2> stats.txt 1> out.fq

# Print reads w/o UMI to stderr
umitools trim --save-invalid unprocessed_fastq NNNNNV 2> invalid.fq 1> out.fq

# Throw error
umitools trim --verbose --save-invalid unprocessed_fastq NNNNNV > out.fq

Both process_fastq and clip_umi would have to be updated to account for this new option. clip_umi would need to return the read instead of just the invalid UMI. If you like this idea, please let me know if you'd like any help implementing and/or testing it.

Paired-end reads?

How does this work with paired-end reads?

Any tutorial?

Hi,

May I ask that is there any tutorial for this tool?

How could I discard reads without pre-specified UMI pattern?

Many thanks!

Best wishes,
Wenhao

not accounting for soft clipping

rmdup is not accounting for soft clipping at each base position. not good.

still not properly accounting for UMIs on the positive strand

append umi info to the left of any spaces in read names

bowtie is clipping read names at spaces, therefore losing umi sequence in the process.

`--mismatches` not working as one could reasonably expect

Given:

ATTG ATTA ATTT AGTA AGTC AGGA

UMI count would be 2: [ATTG, AGTA]

On the other hand, if the order of the UMIs in the file happened to be
arranged differently:

ATTA AGTA AGGA ATTG ATTT AGTC

UMI count would be 3: [ATTA, AGGA, AGTC]

If you only care about unique counts and ordering isn't an issue, then this toolset will work fine. If you want something that works and is much more robust, consider using:

https://github.com/CGATOxford/UMI-tools

Thanks to @jdblischak for pointing this out!

output

Any guidance on interpreting the output?
processing chromosome 1
1 9 10 4 4
1 11 12 2 1
1 29 30 2 2

umitools quantify

It would be useful to have a report / analysis capability in umitools that quantifies issues that arise in UMI analysis.

One question that repeatedly comes up is how much oversampling / saturation of UMIs is happening in a library. For each aligned position, it would be useful to quantify the number of UMIs that are present before and after rmdup.

$ umitools quantify before.rmdup..bam after.rmdup.bam > sample.umi.quant.tab

A potentially useful report of quantification would be a BED3+ file with the umi counts before and after rmdup as the 4th and 5th fields:

chr1 179 180 64 32

A related question is whether 2 libraries prepared with identical UMI tagging strategies can be compared (e.g. per-site) against each other. This would require confirming whether UMIs are saturated at identical, specific positions. This seems potentially difficult and clunky.