rylanci/deduper

Deduper

Part 1

Fork this repo - you should do all your work in your fork of this repository.

Write up a strategy for writing a Reference Based PCR Duplicate Removal tool. That is, given a sam file of uniquely mapped reads, remove all PCR duplicates (retain only a single copy of each read). Develop a strategy that avoids loading everything into memory. You should not write any code for this portion of the assignment. Be sure to:

• Define the problem
• Write examples:
  • Include a properly formated input sam file
  • Include a properly formated expected output sam file
• Develop your algorithm using pseudocode
• Determine high level functions
  • Description
  • Function headers
  • Test examples (for individual functions)
  • Return statement

For this portion of the assignment, you should design your algorithm for single-end data, with 96 UMIs. UMI information will be in the QNAME, like so: NS500451:154:HWKTMBGXX:1:11101:15364:1139:GAACAGGT. Discard any UMIs with errors (or error correct, if you're feeling ambitious).

Part 2

An important part of writing code is reviewing code - both your own and other's. In this portion of the assignment, you will be assigned 3 students' algorithms to review. Be sure to evaluate the following points:

• Does the proposed algorithm make sense to you? Can you follow the logic?
• Does the algorithm do everything it's supposed to do? (see part 1)
• Are proposed functions reasonable? Are they "standalone" pieces of code?

You can find your assigned reviewees on Canvas. You can find your fellow students' repositories at

github.com/<user>/Deduper

Be sure to leave comments on their repositories by creating issues or by commenting on the pull request.

Part 3

Write your deduper function!

Given a SAM file of uniquely mapped reads, remove all PCR duplicates (retain only a single copy of each read). Remember:

• Samtools sort
• Adjust for soft clipping
• Single-end reads
• Known UMIs
• Considerations:
  • Millions of reads – avoid loading everything into memory!
  • Be sure to utilize functions appropriately
  • Appropriately comment code and include doc strings
• CHALLENGE: Include options for
  • Single-end vs paired-end
  • Known UMIs vs randomers
  • Choice of duplicate written to file

You MUST:

• Write Python 3 compatible code
• Include the following argparse options
  • -f, --file: required arg, absolute file path
  • -p, --paired: optional arg, designates file is paired end (not single-end)
  • -u, --umi: optional arg, designates file containing the list of UMIs (unset if randomers instead of UMIs)
  • -h, --help: optional arg, prints a USEFUL help message (see argparse docs)
    • If your script is not capable of dealing with a particular option (ex: no paired-end functionality), your script should print an error message and quit
• Output the first read encountered if duplicates are found
  • You may include an additional argument to designate output of a different read (highest quality or random or ???)
• Output a properly formatted SAM file with “_deduped” appended to the filename
• Name your python script <your_last_name>_deduper.py

deduper's Issues

Peer Review

Hi Ryan,
Awesome work!
Here are my comments for pseudo code:

1. Does the proposed algorithm make sense to you? Can you follow the logic?
   Yes it very well formatted, easy to read and the logic is sensible.

2.Does the algorithm do everything it's supposed to do? (see part 1)
Yes. I think only a minor addition would be to add a proper function header - 'def'

3.Are proposed functions reasonable? Are they "standalone" pieces of code?
Yes. I would not check flag 256, but as it indicates secondary alignment. I think it is 16 which denoted whether the read is forward mapped or reversed.

Again, awesome work!
Pranav

Matt's Peer Review

The algorithm makes sense and is easy to follow. For optimizing the run time, if you extract and check the chromosome first in the pseudocode, I think it removes the need to store it in the condensed_rec list. This is because the list will be flushed for each new one. I was a little confused how “Chroms += Chrom” fit into making the chromosome list, and I didn’t see the list defined elsewhere. Also, the pseudocode shows the UMI being index 2 of qname.split, but I think it is actually index 7.

In parsing the CIGAR string, I think that the forward strand mapped reads should only subtract the soft clipping, and only on the front of the CIGAR string. For reverse, I think that the match/mismatched(M), not included (N), and deletions (D) should be added, as well as soft clipping (S) on the right end of the CIGAR string only. Otherwise I think the pseudocode does everything it is supposed to.

I think it might help the forward or reverse function to return “forward” or “reverse” for the sake of having a good standalone function. Same goes for the mapped/unmapped function. All three functions seem pretty broadly applicable for other software.

Peer_review

Does the proposed algorithm make sense to you? Can you follow the logic?
The propose algorithm makes sense to me and I was able to follow the logic.

Does the algorithm do everything it's supposed to do? (see part 1)
The algorithm seems to do everything that part 1 asked to do.

Are proposed functions reasonable? Are they "standalone" pieces of code?
For your mapped() function, I’m not sure if you need to check for 256. I think that checks for secondary alignment. Also, maybe use regex to address the @ headers at the beginning of a SAM file. I was not aware of the set(), that can definitely be useful.

Overall looks good. Good work.