Covid19 #140
Covid19 #140
Conversation
…undant and adds large amounts of data to batch samples
…undant and adds large amounts of data to batch samples
|
Hi @jrober84, |
bio_hansel/subtyper.py
Outdated
| @@ -228,6 +228,24 @@ def parallel_query_reads(reads: List[Tuple[List[str], str]], | |||
| outputs = [x.get() for x in res] | |||
| return outputs | |||
|
|
|||
| def filter_by_kmer_fraction(df,min_kmer_frac=0.05): | |||
There was a problem hiding this comment.
Hi @jrober84 I'm having trouble understanding what this function is supposed to do.
Is kmer fraction supposed to be like the alternate allele fraction (AF) from variant calling?
What is a noisy kmer? Any kmer that is observed at a low frequency relative to the sum of frequencies of all kmers observed at that position?
kmer_freq / sum(kmer_freqs_at_position) < min_kmer_frac
There was a problem hiding this comment.
Hey @peterk87, I have updated the description in the code to be a bit more clear. But essentially, the function determines the total number of positive and negative kmers covering a position and then determines the percentage of the total coverage each k-mer for the position contributes. In the case where only a single k-mer is present it will always be one and will not ever be filtered. But when both positive and negative are present it will see if the percentage contribution of each k-mer to the total is above the set threshold. If it isn't that specific k-mer gets filtered from the data frame and so the QA/QC module will never see it. I would say this process has some similarity to the AF for variant calling but also will allow the user to configure an acceptable "contamination" level for their sample which for their application is not an issue.
There was a problem hiding this comment.
Hi @jrober84
I'm not sure this code is doing what your description says it's doing. I don't see any accessing of kmer frequency values. If I'm interpreting your description above and in the docstring correctly, you need to be calculating the sum of kmer frequencies at each position. With .value_counts() on refposition you're simply getting a count of how many times you observe each refposition value (similar to Python's Counter). It also doesn't make sense that the refposition count is being divided by the refposition value.
Also, when possible, I think it's a better idea to include all results for the detailed results report for debugging and troubleshooting rather than filtering those results out. I would instead just modify this line:
biohansel/bio_hansel/subtyper.py
Line 289 in d20a00b
Where instead of just getting the subset of results where is_kmer_okay, kmers that pass the kmer fraction threshold could also be filtered for:
df[(df.is_kmer_freq_okay | (df.kmer_fraction >= subtyping_params.min_kmer_fraction))]This way no results are removed from the detailed report, and only the "good" kmers are used to determine the subtype result.
So I would propose adding columns like kmer_fraction and maybe total_refposition_kmer_frequency and does_pass_kmer_fraction_threshold. This would be useful information for potential contamination detection.
There was a problem hiding this comment.
That's a good suggestion, I have updated the code to no longer be filtering the df for failed k-mers. This way all detected k-mers will appear in the detailed report but the QA/QC will be only on the pass list. Pass kmers must pass both the k-mer freq filter and k-mer frac filter. I have added the additional fields that you suggested for troubleshooting purposes.
… being filtered for failed kmers
bio_hansel/subtyper.py
Outdated
| Returns: | ||
| - pd.DataFrame with k-mers with kmer_fraction column | ||
| """ | ||
| position_frequencies = df[['refposition','freq']].groupby(['refposition']).sum().reset_index() |
There was a problem hiding this comment.
It might be easier to use a dict:
position_frequencies = df[['refposition','freq']].groupby(['refposition']).sum().to_dict()The dict should have refposition keys and summed frequency values, so getting total_freq would be easier and clearer:
total_freq = position_frequencies[row.refposition]
bio_hansel/subtyper.py
Outdated
| position_frequencies = df[['refposition','freq']].groupby(['refposition']).sum().reset_index() | ||
| percentages = [] | ||
| total_refposition_kmer_frequencies = [] | ||
| for index,row in df.iterrows(): |
There was a problem hiding this comment.
I'd recommend using .itertuples() for performance reasons:
for row in df.itertuples():
refposition = row.refposition # cannot do string based access (row['refposition']) with tuplesYou could also look into using apply instead of using a for-loop:
def get_kmer_fraction(row):
total_freq = position_frequencies.get(row.refposition, 0)
return row.freq / total_freq if total_freq > 0 else 0.0
df['kmer_fraction'] = df.apply(get_kmer_fraction, axis=1)
df['total_refposition_kmer_frequency'] = df.apply(lambda row: position_frequencies.get(row.refposition, 0), axis=1)
There was a problem hiding this comment.
Hi @jrober84
Thanks for making the suggested changes and addressing my comments and questions.
The results I'm getting from Nanopore and Illumina data make sense with the low abundance kmer matches excluded from subtype calling yet still present in the detailed report. So everything looks good on my end and ready to merge into dev.
This update is to make some minor modifications to BioHansel to be more compatible with amplicon based SARS-COV-2 data by introducing the min-kmer-frac parameter which works in combination with the min-kmer-cov parameter to ignore k-mers present in a sample which are less than the defined percentage set by the program with a default of 0.05. Additionally, the max-kmer-cov has been increased to not cause conflicts with high covered regions being ignored. The kmer report has also removed the error message field since this adds a lot of extra repeated data for these reports which adds up for larger sample pools.