<commit_before>#!/usr/bin/env python
# -*- coding: utf-8 -*-

# ------------------------------------------------------------------------------
# 
# Author: <NAME>
# Email: <EMAIL>
# Version: 1.2.0
# Date: 20170628
# Project: RNA FISH oligo design
# Description:	filter BLASTN output based on:
# 					- homology (#PM/k),
# 					- number of OT, and
# 					- common (saturated) OTs.
# 
# Notes:
# 		PM: Perfect Match.
# 		OT: Off Target.
# 		The script should be run on single-gene outputs.
# 
# Changelog:
# 		1.0.0: first implementation.
# 		1.1.0: added argparser support.
# 		1.2.0: changed saturation OT filter.
# 
# ------------------------------------------------------------------------------



# DEPENDENCIES =================================================================

import argparse
import numpy as np

# PARAMETERS ===================================================================


# Add script description
parser = argparse.ArgumentParser(
	description = 'Filter BLASTN output.'
)

# Add mandatory arguments
parser.add_argument('blastInput', type = str, nargs = 1,
	help = 'Path to BLAST input fasta file.')
parser.add_argument('blastOutput', type = str, nargs = 1,
	help = 'Path to BLAST output with outfmt 6.')
parser.add_argument('geneTranscriptTable', type = str, nargs = 1,
	help = """Path to table with TRANSCRIPT_ID:GENE_SYMBOL
	tabulation-separated columns.""", default = [30])
parser.add_argument('output', type = str, nargs = 1,
	help = 'Path to output fasta file, after filtering.')

# Add arguments with default value
parser.add_argument('-k', type = int, nargs = 1,
	metavar = 'k', help = """Oligonucleotide length in nt.
	Default: 30.""", default = [30])
parser.add_argument('-t', '--homology-thr', type = float, nargs = 1,
	metavar = 'ht', help = """Threshold on maximum homology, as fraction of k.
	Accepts float values from 0 to 1.
	Default: .85""", default = [.85])
parser.add_argument('-e', '--gene-thr', type = int, nargs = 1,
	metavar = 'gt', help = """Threshold on the number of off-targets gene,
	for a single oligo.
	Default: 20""", default = [20])
parser.add_argument('-s', '--saturation-thr', type = int, nargs = 1,
	metavar = 'st', help = """
	Threshold on the number of oligos off-targeting a gene,
	for the selections of 'saturated' off-target genes. Oligos targeting
	a saturated off-target are filtered out.
	Default: 5""", default = [5])

# Parse arguments
args = parser.parse_args()

# Assign to in-script variables
blast_input = args.blastInput[0]
blast_output = args.blastOutput[0]
gene_transcript_table = args.geneTranscriptTable[0]
output_file = args.output[0]
k = args.k[0]
homThr = args.homology_thr[0]
gene_ot_thr = args.gene_thr[0]
oligo_ot_thr = args.saturation_thr[0]

# Log to screen the settings
print("""
Settings:
              BLASTN input : %s
             BLASTN output : %s
     Gene-Transcript table : %s
               Output file : %s
                         K : %d
        Homology threshold : %f
       Gene Off-Target thr : %d
 Off-Target saturation thr : %d

""" % (blast_input, blast_output, gene_transcript_table, output_file,
	k, homThr, gene_ot_thr, oligo_ot_thr))


# FUNCTIONS ====================================================================

# RUN ==========================================================================
print("Run:")

# Build TRANSCRIPT_ID:GENE_SYMBOL dictionary -----------------------------------
print(" · Building TRANSCRIPT_ID:GENE_SYMBOL dictionary...")

# Initialize empty TRANSCRIPT_ID:GENE_SYMBOL dictionary
trn_gene_dict = {}

# Read table line by line
with open(gene_transcript_table) as gttf:
	for line in gttf:

		# Split every line in two fields
		tmp = line.strip().split('\t')

		# Store key,value couples in the dictionary
		trn_gene_dict[tmp[0]] = tmp[1]

# Work on BLASTN output --------------------------------------------------------

# Filter based on max homology
# ----------------------------
print(" · Filtering based on maximum homology...")

# Initalize empty maximum homology dictionary
max_homology = {}

# Read table line by line
with open(blast_output) as bof:
	for line in bof:

		# Split every line by column
		tmp = line.strip().split('\t')

		# Identify oligomer ID
		OID = tmp[0].split(':')[0]

		# Identify target gene
		target = tmp[0].split('_')[0]

		# Identify transcript ID
		transcript = tmp[1].split('.')[0]

		# Calculate homology
		homology = (int(tmp[3]) - int(tmp[4])) / float(k)

		# If it is an off-target
		if target != trn_gene_dict[transcript]:
			if not OID in max_homology.keys():
				max_homology[OID] = homology
			else:
				# With homology higher than the threshold
				if homology > max_homology[OID]:
					max_homology[OID] = homology

# Identify oligos that pass the threshold
pass_homology = []
for OID in max_homology.keys():
	if max_homology[OID] < homThr:
		pass_homology.append(OID)
# Log
print(" >>> %d oligos do not have any off-targets." % (len(pass_homology),))
print(" >>> %d oligos have off-targets." % (len(max_homology),))
print(" >>> Saving off-target free oligos. Analyzing further the rest.")
<commit_msg>Implemented default max homology of 0 for oligos with no off-targets. Fixed log.<commit_after>#!/usr/bin/env python
# -*- coding: utf-8 -*-

# ------------------------------------------------------------------------------
# 
# Author: <NAME>
# Email: <EMAIL>
# Version: 1.2.0
# Date: 20170628
# Project: RNA FISH oligo design
# Description:	filter BLASTN output based on:
# 					- homology (#PM/k),
# 					- number of OT, and
# 					- common (saturated) OTs.
# 
# Notes:
# 		PM: Perfect Match.
# 		OT: Off Target.
# 		The script should be run on single-gene outputs.
# 
# Changelog:
# 		1.0.0: first implementation.
# 		1.1.0: added argparser support.
# 		1.2.0: changed saturation OT filter.
# 
# ------------------------------------------------------------------------------



# DEPENDENCIES =================================================================

import argparse
import numpy as np

# PARAMETERS ===================================================================


# Add script description
parser = argparse.ArgumentParser(
	description = 'Filter BLASTN output.'
)

# Add mandatory arguments
parser.add_argument('blastInput', type = str, nargs = 1,
	help = 'Path to BLAST input fasta file.')
parser.add_argument('blastOutput', type = str, nargs = 1,
	help = 'Path to BLAST output with outfmt 6.')
parser.add_argument('geneTranscriptTable', type = str, nargs = 1,
	help = """Path to table with TRANSCRIPT_ID:GENE_SYMBOL
	tabulation-separated columns.""", default = [30])
parser.add_argument('output', type = str, nargs = 1,
	help = 'Path to output fasta file, after filtering.')

# Add arguments with default value
parser.add_argument('-k', type = int, nargs = 1,
	metavar = 'k', help = """Oligonucleotide length in nt.
	Default: 30.""", default = [30])
parser.add_argument('-t', '--homology-thr', type = float, nargs = 1,
	metavar = 'ht', help = """Threshold on maximum homology, as fraction of k.
	Accepts float values from 0 to 1.
	Default: .85""", default = [.85])
parser.add_argument('-e', '--gene-thr', type = int, nargs = 1,
	metavar = 'gt', help = """Threshold on the number of off-targets gene,
	for a single oligo.
	Default: 20""", default = [20])
parser.add_argument('-s', '--saturation-thr', type = int, nargs = 1,
	metavar = 'st', help = """
	Threshold on the number of oligos off-targeting a gene,
	for the selections of 'saturated' off-target genes. Oligos targeting
	a saturated off-target are filtered out.
	Default: 5""", default = [5])

# Parse arguments
args = parser.parse_args()

# Assign to in-script variables
blast_input = args.blastInput[0]
blast_output = args.blastOutput[0]
gene_transcript_table = args.geneTranscriptTable[0]
output_file = args.output[0]
k = args.k[0]
homThr = args.homology_thr[0]
gene_ot_thr = args.gene_thr[0]
oligo_ot_thr = args.saturation_thr[0]

# Log to screen the settings
print("""
Settings:
              BLASTN input : %s
             BLASTN output : %s
     Gene-Transcript table : %s
               Output file : %s
                         K : %d
        Homology threshold : %f
       Gene Off-Target thr : %d
 Off-Target saturation thr : %d

""" % (blast_input, blast_output, gene_transcript_table, output_file,
	k, homThr, gene_ot_thr, oligo_ot_thr))


# FUNCTIONS ====================================================================

# RUN ==========================================================================
print("Run:")

# Build TRANSCRIPT_ID:GENE_SYMBOL dictionary -----------------------------------
print(" · Building TRANSCRIPT_ID:GENE_SYMBOL dictionary...")

# Initialize empty TRANSCRIPT_ID:GENE_SYMBOL dictionary
trn_gene_dict = {}

# Read table line by line
with open(gene_transcript_table) as gttf:
	for line in gttf:

		# Split every line in two fields
		tmp = line.strip().split('\t')

		# Store key,value couples in the dictionary
		trn_gene_dict[tmp[0]] = tmp[1]

# Work on BLASTN output --------------------------------------------------------

# Filter based on max homology
# ----------------------------
print(" · Filtering based on maximum homology...")

# Initalize empty maximum homology dictionary
max_homology = {}

# Read table line by line
with open(blast_output) as bof:
	for line in bof:

		# Split every line by column
		tmp = line.strip().split('\t')

		# Identify oligomer ID
		OID = tmp[0].split(':')[0]

		# Identify target gene
		target = tmp[0].split('_')[0]

		# Identify transcript ID
		transcript = tmp[1].split('.')[0]

		# Calculate homology
		homology = (int(tmp[3]) - int(tmp[4])) / float(k)

		# If it is an off-target
		if target != trn_gene_dict[transcript]:
			if not OID in max_homology.keys():
				max_homology[OID] = homology
			else:
				# With homology higher than the threshold
				if homology > max_homology[OID]:
					max_homology[OID] = homology
		elif not OID in max_homology.keys():
			max_homology[OID] = 0

# Identify oligos that pass the threshold
pass_homology = []
for OID in max_homology.keys():
	if max_homology[OID] < homThr:
		pass_homology.append(OID)

# Log
print(" >>> %d oligos do not have any off-targets." % (len(pass_homology),))
print(" >>> %d oligos have off-targets." % (len(max_homology)-len(pass_homology),))
print(" >>> Saving off-target free oligos. Analyzing further the rest.")

{'fraction_non_alphanumeric': 0.22313569977754136, 'fraction_numerical': 0.01296063449076313, 'mean_word_length': 3.0901898734177213, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6617810', 'n_tokens_mistral': 3587, 'n_tokens_neox': 3180, 'n_words': 983}