From 3685d4a760bdf3181588c523c4ef5c278a35d3b5 Mon Sep 17 00:00:00 2001 From: Runxi Shen Date: Mon, 12 Apr 2021 13:20:23 +0800 Subject: [PATCH] update celescope/tools/consensus.py. Incorporate sequencing error into consensus identification. --- celescope/tools/consensus.py | 105 ++++++++++++++++++++++++----------- 1 file changed, 72 insertions(+), 33 deletions(-) diff --git a/celescope/tools/consensus.py b/celescope/tools/consensus.py index 3ada98fb..114b1654 100644 --- a/celescope/tools/consensus.py +++ b/celescope/tools/consensus.py @@ -3,14 +3,22 @@ import gzip import numpy as np import subprocess import os +import pickle +import scipy.stats as sp_stats from xopen import xopen -from collections import defaultdict +from collections import defaultdict, OrderedDict, Counter from itertools import groupby +from scipy.special import softmax from celescope.tools.utils import * from celescope.tools.report import reporter from celescope.tools.Reporter import Reporter +### CONSTANTS +## prob threshold for calling a nucleotide as ambiguous (detailed calculation in codes below) +AMBIGUOUS_NUC_CONFIDENCE_THRES = 0.9 + + @add_log def sort_fastq(fq, fq_tmp_file, outdir): tmp_dir = f'{outdir}/tmp' @@ -22,16 +30,18 @@ def sort_fastq(fq, fq_tmp_file, outdir): @add_log -def sorted_dumb_consensus(fq, outfile, threshold): +def sorted_dumb_consensus(fq, outfile, umi_reads_count_pckl, threshold): ''' consensus read in name-sorted fastq - output (barcode,umi) consensus fastq + store the umi counts in a dict + output (barcode,umi) consensus fastq and umi-count dict ''' read_list = [] n_umi = 0 total_ambiguous_base_n = 0 length_list = [] out_h = xopen(outfile, 'w') + umi_reads_count_dict = OrderedDict() def keyfunc(read): attr = read.name.split('_') @@ -48,11 +58,15 @@ def sorted_dumb_consensus(fq, outfile, threshold): prefix = "_".join([barcode, umi]) read_name = f'{prefix}_{n_umi}' out_h.write(fastq_line(read_name, consensus_seq, consensus_qual)) + umi_reads_count_dict[(barcode, umi)] = len(read_list) if n_umi % 10000 == 0: sorted_dumb_consensus.logger.info(f'{n_umi} UMI done.') total_ambiguous_base_n += ambiguous_base_n length_list.append(con_len) + with open(umi_reads_count_pckl, 'wb') as handle: + pickle.dump(umi_reads_count_dict, handle, protocol=pickle.HIGHEST_PROTOCOL) + out_h.close() return n_umi, total_ambiguous_base_n, length_list @@ -62,30 +76,42 @@ def wrap_consensus(fq, outdir, sample, threshold): fq_tmp_file = f'{outdir}/{sample}_sorted.fq.tmp' sort_fastq(fq, fq_tmp_file, outdir) outfile = f'{outdir}/{sample}_consensus.fq' + umi_reads_count_pckl = f'{outdir}/{sample}_umi_read_counts.pckl' n, total_ambiguous_base_n, length_list = sorted_dumb_consensus( - fq=fq_tmp_file, outfile=outfile, threshold=threshold) + fq=fq_tmp_file, outfile=outfile, + umi_reads_count_pckl=umi_reads_count_pckl, + threshold=threshold) return outfile, n, total_ambiguous_base_n, length_list -def dumb_consensus(read_list, threshold=0.5, ambiguous='N', default_qual='F'): +def dumb_consensus(read_list, threshold=0.5, ambiguous='N'): ''' - This is similar to biopython dumb_consensus. It will just go through the sequence residue by residue and count up the number of each type - of residue (ie. A or G or T or C for DNA) in all sequences in the - alignment. If the percentage of the most common residue type is + of residue (ie. A or G or T or C for DNA) in all reads in the alignment. + For each position in each read, we will calculate a probability of sequencing confidence by + using the number of occurrences of the nucleotide (k) and the PHRED score of the nucleotide + representing the error probability (phred) of sequencing by k*(1-phred), and calculate the + probability of observing the nucleotide at the position by normalizing the counts to an array + of probabilities summing to 1. The final consensus would be the nuc with the largest prob. + + If the probability of observing the most likely residue type is greater then the passed threshold, + then we will add that residue type and its newly calculated quality score converted to ASCII code + using the same procedure as PHRED chr(Q=-10*log10(1-P)+33). If the percentage of the most common residue type is greater then the passed threshold, then we will add that residue type, otherwise an ambiguous character will be added. elements of read_list: [entry.sequence,entry.quality] ''' con_len = get_read_length(read_list, threshold=threshold) - consensus_seq = "" + consensus = "" consensus_qual = "" ambiguous_base_n = 0 + for n in range(con_len): - atom_dict = defaultdict(int) - quality_dict = defaultdict(int) - num_atoms = 0 + total_obs = 0 + count_atom_occur = defaultdict(int) + atom_phred_error = defaultdict(float) + atom_qual = defaultdict(str) for read in read_list: # make sure we haven't run past the end of any sequences # if they are of different lengths @@ -93,30 +119,43 @@ def dumb_consensus(read_list, threshold=0.5, ambiguous='N', default_qual='F'): quality = read[1] if n < len(sequence): atom = sequence[n] - atom_dict[atom] += 1 - num_atoms = num_atoms + 1 - base_qual = quality[n] - quality_dict[base_qual] += 1 - + ## convert the phred qual to error prob + phred_2_error = 10**-((ord(base_qual)-33)/10.0) + count_atom_occur[atom] += 1 + if (phred_2_error= num_atoms * threshold: - consensus_atom = atom - break - if consensus_atom == ambiguous: + ## calculate prob of obs using occurrence*(1-phred) per nuc + ## the phred error is most accurate call of this nucleotide + if (len(atom_phred_error.items()) == 1): + atom, phred_error = list(atom_phred_error.items())[0] + max_prob = 1 - phred_error + else: + atom_phred_list = list(atom_phred_error.items()) + atom_occurrences = [count_atom_occur[atom]*(1-phred) for atom, phred in atom_phred_list] + softmax_prob = softmax(atom_occurrences) + max_prob = softmax_prob[np.argmax(softmax_prob)] + atom = atom_phred_list[np.argmax(softmax_prob)][0] + + ## check whether prob passes AMBIGUOUS_NUC_CONFIDENCE_THRES + if (max_prob > AMBIGUOUS_NUC_CONFIDENCE_THRES): + consensus_atom = atom + else: ambiguous_base_n += 1 - consensus_seq += consensus_atom - - max_freq_qual = 0 - consensus_base_qual = default_qual - for base_qual in quality_dict: - if quality_dict[base_qual] > max_freq_qual: - max_freq_qual = quality_dict[base_qual] - consensus_base_qual = base_qual - - consensus_qual += consensus_base_qual - return consensus_seq, consensus_qual, ambiguous_base_n, con_len + consensus += consensus_atom + + ## if the max prob is larger than the best possibility + ## we keep the best possibility among all reads + if (max_prob >= 1-atom_phred_error[atom]): + consensus_qual += atom_qual[atom] + else: + consensus_qual += chr(int(-10*np.log10(1-max_prob)+33)) + + return consensus, consensus_qual, ambiguous_base_n, con_len def get_read_length(read_list, threshold=0.5): -- Gitee