Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Working #1

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Working #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b64cfa1
Jupyter Notebook Added to GitHub.com
MaxHills Oct 23, 2018
07d8428
Add test SAM files: INPUT/OUTPUT
MaxHills Oct 23, 2018
a13aa6a
test file for running jupyter notebook
MaxHills Oct 23, 2018
42b971d
Update final version of jupyter notebook w/ corrections
MaxHills Oct 23, 2018
d116c3b
"One last update..."
MaxHills Oct 23, 2018
922f246
Add Hills_deduper.py
MaxHills Nov 6, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
262 changes: 262 additions & 0 deletions .ipynb_checkpoints/Deduper - Part 1-checkpoint.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,262 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Deduper \n",
"--- \n",
"## Part 1 \n",
"--- \n",
"### Define the Problem \n",
"--- \n",
"- _PCR duplicates arise during PCR amplification of original fragmented, adapter-ligated DNA. They become problematic when 2 copies of the same molecule get onto different primer lawns on an Illumina flowcell. PCR-duplicate removal (or tagging) is necessary to mitigate PCR bias, which is introduced when certain molecules ligate adapters more favorably and are disproportionately amplified. The greater numbers of these preferentially amplified molecules makes them more likely to contaminate neighboring lawns on the flowcell. Failure to remove these duplicates could limit your ability to draw conclusions from your data, as biological information may be confounded by PCR bias._ \n",
"\n",
"### Write Examples \n",
"--- \n",
" - Input SAM file: \n",
" - **sample\\_input.sam** \n",
" - 6 conditions for forward read:\n",
" - Initial UMI (non-duplicate)\n",
" - Same UMI, different chromosome\n",
" - Same UMI, different POS \n",
" - Same UMI, different strand\n",
" - Same UMI, DUPLICATE w/ soft-clipping\n",
" - Ssme UMI, DUPLICATE w/o soft-clipping\n",
" - 6 conditions for reverse read:\n",
" - Same as above\n",
" - Output SAM file: \n",
" - **sample\\_output.sam** \n",
" - 8 reads: 4 from forward, 4 from reverse.\n",
"\n",
"### Develop Algorithm w/ Pseudocode \n",
"--- "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#################\n",
"## PSEUDO-CODE ##\n",
"#################\n",
"# To run SAMtools on the command line from a Python script...\n",
"from subprocess import call \n",
"# Import regular expressions.\n",
"import re\n",
"\n",
"# run SAMtools sort on file(s).\n",
"call(\"samtools sort -l 0 -m maxMem -o <out.sam> -T <out.prefix> <input_file.sam>\", shell=True)\n",
"\n",
"# open the file produced from 'samtools sort' (above) as fin, open the output file for THIS program as fout.\n",
"# open the file of known barcodes as index.\n",
"with open(\"out.sam\") as fin, with open(\"dedup_output.sam\",'a') as fout, with open(\"STL96.txt\") as index:\n",
" # Create an empty set to hold real indexes.\n",
" valid_IDs = set()\n",
" # Add each index to valid_IDs.\n",
" for i in index:\n",
" i = i.strip()\n",
" valid_IDs.add(i)\n",
" # Create a counter for PCR-duplicates.\n",
" PCR_dup = 0\n",
" # Create a dictionary to keep track of UMI information.\n",
" UMI_dict = {}\n",
" # Keep track of total reads.\n",
" tot_rds = 0\n",
" # Iterate through each line of the input SAM file.\n",
" for line in fin:\n",
" # If the line is a header line...\n",
" if re.match(r'^@[A-Z]{2}\\s', line):\n",
" # Store the line w/newline in the temporary string 'header'.\n",
" header = (line + \"\\n\")\n",
" # Write the header to fout.\n",
" fout.write(header)\n",
" # Continue iterating thru 'fin'.\n",
" continue\n",
" else:\n",
" # Increment the number of total reads.\n",
" tot_rds += 1\n",
" # Get the relevant information to check for PCR duplicates.\n",
" parse_read = parseSAM(line)\n",
" # Store info in separate variables.\n",
" this_UMI = parse_read[0]\n",
" this_chr = parse_read[1]\n",
" this_pos = parse_read[2]\n",
" this_CIGAR = parse_read[3]\n",
" this_strand = parse_read[4]\n",
" # If current UMI is in valid_IDs...\n",
" if this_UMI in valid_IDs:\n",
" # If current read HASN'T undergone soft-clipping...\n",
" if this_pos == truePOS(this_pos,this_CIGAR):\n",
" # And current read is not in the dictionary.\n",
" if this_UMI not in UMI_dict:\n",
" # Add it to the dictionary.\n",
" UMI_dict[this_UMI] = (this_chr,this_pos,this_strand)\n",
" read = line + \"\\n\"\n",
" # Write the line to fout.\n",
" fout.write(read)\n",
" # If the UMI is already in the dictionary, check for PCR duplicate...\n",
" elif:\n",
" if this_chr == UMI_dict[this_UMI][0] and this_pos == UMI_dict[this_UMI][1] and this_strand == UMI_dict[this_UMI[2]]:\n",
" # Same UMI, same chromosome, same position, and same strand... This is a PCR duplicate.\n",
" # Don't write the line to file. Continue through file.\n",
" # Increment PCR-duplicate counter.\n",
" PCR_dup += 1\n",
" continue\n",
" # If the UMI is in dictionary, but not a duplicate...\n",
" else:\n",
" read = line + \"\\n\"\n",
" # Write the read to fout.\n",
" fout.write(read)\n",
" # The current read HAS undergone soft-clipping...\n",
" else:\n",
" # Correct the position.\n",
" this_pos = truePOS(this_pos,this_CIGAR)\n",
" # Check for PCR-duplicate.\n",
" if this_chr == UMI_dict[this_UMI][0] and this_pos == UMI_dict[this_UMI][1] and this_strand == UMI_dict[this_UMI[2]]:\n",
" # Same UMI, same chromosome, same position, and same strand... This is a PCR duplicate.\n",
" # Don't write the line to file. Continue through file.\n",
" # Increment PCR-duplicate counter.\n",
" PCR_dup += 1\n",
" continue\n",
" # The read is not a PCR duplicate, write to fout.\n",
" else:\n",
" read = line + \"\\n\"\n",
" fout.write(read)\n",
" # If the UMI isn't valid, DO NOT write the current read to the output file. \n",
" else:\n",
" # Error correction on UMIs would occur in this block.\n",
" continue\n",
" \n",
" \n",
" \n",
" \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Determine High-level functions \n",
"--- "
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"import re\n",
"\n",
"def truePOS(POS,CIGAR):\n",
" \"\"\"\n",
" description: truePOS checks if soft-clipping has occurred at the left-most mapped position on a read.\n",
" If soft-clipping has occurred, truePOS performs a calculation to determine the INT value for the original \n",
" left-most-mapped position and returns this corrected pos. If soft-clipping has not occurred, truePOS returns\n",
" the original INT value of POS.\n",
" INPUT: POS = A STR representing the left-most-mapped position of the read.\n",
" CIGAR = A STR representing the CIGAR-string for the read.\n",
" OUTPUT: pos = A SC-corrected value for POS in reads where soft-clipping has occurred.\n",
" POS = The original integer value for POS\n",
" \"\"\"\n",
" # Make POS an INT for calculations.\n",
" POS = int(POS)\n",
" # Search for left-most soft-clipping.\n",
" sc = re.search(\"([0-9]{1,2})S[0-9]{1,2}\", CIGAR)\n",
" # If it's found...\n",
" if sc:\n",
" # Set the difference to the INT value before the first S, for calculations.\n",
" diff = int(sc.group(1))\n",
" # Create the new corrected position, to check for PCR duplicates.\n",
" pos = POS - diff\n",
" # Return the corrected position.\n",
" return pos\n",
" # If there's no soft-clipping...\n",
" else:\n",
" # Return the original position\n",
" return POS\n",
" \n",
"def parseSAM(lineSAM):\n",
" \"\"\"\n",
" description: parseSAM takes 1 line from a SAM file and returns the UMI, chromosome #, POS, CIGAR, and the strand.\n",
" INPUT: A SAM file read-line\n",
" OUTPUT: UMI,CHROM,POS,CIGAR,STRAND\n",
" \"\"\"\n",
" umi = re.search(\".*:([A,C,G,T]{8})\\t\", lineSAM)\n",
" barcode = umi.group(1)\n",
" line = lineSAM.split(\"\\t\")\n",
" chromo = line[2]\n",
" POS = line[3]\n",
" CIGAR = line[5]\n",
" flag = int(line[1])\n",
" if((flag & 16) != 16):\n",
" strand = \"+\"\n",
" else:\n",
" strand = \"-\"\n",
" return barcode,chromo,POS,CIGAR,strand"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"# Test truePOS()\n",
"assert(truePOS(\"500\",\"10S90M1S\")) == 490\n",
"assert(truePOS(\"500\",\"90M11S\")) == 500\n",
"assert(truePOS(\"500\",\"1S98M2S\")) == 499\n",
"assert(truePOS(\"1\",\"101M\")) == 1"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"# Test parseSAM()\n",
"import re\n",
"# Use a known test file on local computer to test parseSAM().\n",
"with open(\"testset3.sam\") as f:\n",
" header = \"\"\n",
" LN = 0\n",
" for line in f:\n",
" if re.match(r'^@[A-Z]{2}\\s', line):\n",
" continue\n",
" else:\n",
" LN += 1\n",
" if LN < 9:\n",
" stuff = parseSAM(line)\n",
" assert(len(stuff[0])) == 8\n",
" assert(stuff[1]) == \"1\"\n",
" assert(stuff[3]) == \"71M\"\n",
" assert(stuff[4]) == \"+\""
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Loading