Jessime
commented
3 years ago Hey Mauro, I've got lots of thoughts about this one!
At the highest level, I think you're reaching the limit of what it makes sense to expose at the command line, as opposed to writing some simple python scripts and having the flexibility that comes along with that.
by specifying custom alphabets
This is neat! Did you know that I've already got support for custom alphabets? You can see on this line, that "AGTC" is just the default. You could count protein kmers if you wanted, haha.
Inside of a python script, this would just be, for example
counter = BasicCounter('gencode.fa', alphabet='AGTCWS')
counter.get_counts()input flag that would cause them to exclude k-mers with "0" counts across all counted RNAs.
The generalization of what you're asking here is, "Can we filter kmers?"
As I'm sure you remember, a big surprise for us was that we couldn't find a subset of kmers that gave us more predictive power than the full set of 6mers. Part of the reason it was so surprising was because I spent a ton of time trying. You've listed a few here:
max(count_column) == 0max(count_column) < thresholdvariance(count_column) < threshold
I can think of several dozen more that we tried over the years. A few were:
- Picking at random
- GC content of kmer
- Covariance with other kmers
- Mean of counts
- More complicated machine learning attempts
I say all of this to say, I explicitly chose not to build something like seekr_filter_counts because there are so many options and it wasn't clear which of them were worth the effort to build into the official package.
On the other hand, a single filter is pretty easy from within python. if you've run seekr_count_kmers, you could do:
def has_vals(col):
return max(col) > 0 # your custom filter
df = pd.read_csv('kmers.csv', index_col=0) # read
filtered_df = df.loc[df.apply(has_vals)] # apply filter
filtered_df.to_csv('filtered_kmers.csv') # writewhich would exclude k-mers with zero counts. And notice that only the line commented as # your custom filter is actually interesting.
What might make sense is to create a seekr_filter_counts that allows users to specify an arbitrary filter function. That'd allow you to incorporate filtering into your your bash scripts, but wouldn't make it necessary to curate all of your favorite filters and go through me each time you want to try a new filter. This would still require someone to write the python filters, of course, but it sounds like the ones you're thinking about right now are easy enough that someone in lab could write them up without too much trouble.
To summarize:
- There are 100s of filters one might apply; it probably doesn't make sense to start adding them to
seekrright now. - It's easy to write one off scripts to filter columns.
- If you're really interested in having filtering as part of the bash scripting pipeline, we should talk more.
Hope this is helpful! 😄
jmcalabr
dan-sprague
This is not really an issue, just a suggestion -- hope all is well!
We are finding new ways to use SEEKR by specifying custom alphabets. In these experiments, we are encountering words or k-mers that never occur, even in the entire mouse or human transcriptome.
Here (and probably in other use cases as well), it would be useful if seekr_norm_vectors and seekr_count_kmers could take an input flag that would cause them to exclude k-mers with "0" counts across all counted RNAs. Or even a flag that would allow exclusion of k-mers that do not exceed a threshold of variance across the counted sets of RNAs.
Thank you and take care, Mauro