Selene is a Python library and command line interface for training deep neural networks from biological sequence data such as genomes.

Added method and unit tests for creating in silico mutagenesis result matrix for one feature as described in #27

This matrix is then passed to seaborn and plotted as a heatmap.

Hi,

I want to test the trained model in https://github.com/FunctionLab/selene/blob/master/tutorials/quickstart_training/quickstart_training.ipynb, the used .yml file is shown below:

ops: [evaluate] model: { # TODO: update this line with the absolute path to the file. path: /root/selene/learn/deeperdeepsea.py, class: DeeperDeepSEA, class_args: { sequence_length: 1000, n_targets: 1, }, non_strand_specific: mean }

sampler: !obj:selene_sdk.samplers.file_samplers.BedFileSampler { filepath: /root/selene/learn/training_outputs/test_data.bed, reference_sequence: !obj:selene_sdk.sequences.Genome { input_path: /root/selene/learn/male.hg19.fasta }, n_samples: 1200, sequence_length: 1000, targets_avail: True, n_features: 2, }

evaluate_model: !obj:selene_sdk.EvaluateModel { features: !obj:selene_sdk.utils.load_features_list { input_path: /root/selene/learn/distinct_features.txt }, trained_model_path: /root/selene/learn/training_outputs/best_model.pth.tar, n_test_samples: 1200, batch_size: 64, report_gt_feature_n_positives: 10, use_cuda: True,
}

random_seed: 1447 output_dir: /root/selene/learn/training_outputs create_subdirectory: False

Then I got this error:

Traceback (most recent call last): File "main.py", line 12, in parse_configs_and_run(configs) File "/root/anaconda3/lib/python3.7/site-packages/selene_sdk/utils/config_utils.py", line 341, in parse_configs_and_run execute(operations, configs, current_run_output_dir) File "/root/anaconda3/lib/python3.7/site-packages/selene_sdk/utils/config_utils.py", line 208, in execute evaluate_model = instantiate(evaluate_model_info) File "/root/anaconda3/lib/python3.7/site-packages/selene_sdk/utils/config.py", line 239, in instantiate return _instantiate_proxy_tuple(proxy, bindings) File "/root/anaconda3/lib/python3.7/site-packages/selene_sdk/utils/config.py", line 144, in _instantiate_proxy_tuple obj = proxy.callable(**kwargs) File "/root/anaconda3/lib/python3.7/site-packages/selene_sdk/evaluate_model.py", line 134, in init if type(self.reference_sequence) == Genome and
AttributeError: 'EvaluateModel' object has no attribute 'reference_sequence'

Would you like to help me debug?

Thank you very much!

Yours sincerely, Shusen

Hi, I am new to your tool and after going through your case studies I got a general overview of the framework. I just want to ask can I supply peaks from different histone modifications ChIP-Seq (like H3K27ac, H3K9ac) as features for a particular tissue and then use them to train the model?

Hi, I have been reading the documentation and I'm still not sure what the output scores for getting predictions from a trained model means. I noticed that the scores are all from 0-1, is it the probability that a TF will bind to an input region and what is this probability based on?

Another question I have is that if I set my "center_bin_to_predict" to be 200 when training the model, and my "feature_thresholds" is 0.5, do my input TF binding regions have to be at least 200bp long for Selene to classify it as a "binding region"

thanks!

Michelle

Currently, I think this issue is only applicable to the sampler module. We may eventually extend the same idea to other modules in selene.

A user should be able to choose what kind of data sampler they want to use when training their model. Each of these samplers has very different input requirements. For example, an IntervalsSampler specific config file would have a lot more parameters compared to a MatFilesSampler config file. In the YAML file, a user should specify the type of sampler they want to use. The value of type should be the class name of the sampler. The remaining parameters could be specified under a parameter key:

type: IntervalsSampler
parameters:
    genome: <path to indexed genome fasta file>
    query_feature_data: <path to tabix-indexed features .bed.gz file>
    ...

Thanks @bmacedo-lgtm for catching this issue:

it turns out that the YAML parser we are using will only accept a single custom metrics function with the present syntax described in the Selene CLI docs.

metrics: {
     roc_auc: !import:sklearn.metrics.roc_auc_score,
     average_precision: !import:sklearn.metrics.average_precision_score
 },

that is, the above snippet that is described in our docs will crash because the parser expects every !import:<method/module> tag to come with a dictionary of input arguments, even if the dictionary is empty.

It seems like Selene will run with no errors if I just modify this to

metrics: {
     roc_auc: !import:sklearn.metrics.roc_auc_score {},
     average_precision: !import:sklearn.metrics.average_precision_score {}
 },

but more testing is needed to better characterize and handle this bug.

I followed the getting started tutorial and trained the GM12878 CTCF binding model. Then I inputted enhancer regions and compared the model's predictions to experimental binding of CTCF in these regions. The model's predictions were vastly different from the experimental CHIP seq binding and I'm wondering why that is?

best,

Michelle

Hi,

Thanks for developing this package. Its very useful for people who wants to get started with machine learning in genomics.

I have a very naive question about Getting Started tutorial. I would like to know whats the specific role of "deepsea_TF_intervals.txt" file in training the model. Are you using the deepsea_TF_intervals file to predict the CTCF binding sites ? or the model is trained on certain number of CTCF peaks and tries to predict the rest ? in this case whats the role of deepsea_TF_intervals ?

Thanks, Goutham A

PS: Is there any google groups to ask some methodological or theoretical questions ?

We would like to include a module or function that allows users to visualize the effects of single base substitutions. This would be a heatmap of the changes in model prediction for a single genomic feature, where the columns of the heatmap correspond to each base position in the input sequence. The rows of the heatmap are the base substitutions.

In in silico mutagenesis, the rows would be all 4 bases (with the reference/original base grayed out) because we substitute in all other bases.

@jzthree I forgot what we decided for visualization of variant effect prediction. Would you be able to provide a description here?

The recent addition of

torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False

has led to a 10x slowdown of variant effect prediction (e.g. taking 1000s to process 10k variants when previously it took 100s to process 10k variants) for one of my DL models. I've tested that this is the reason for the slowdown by running the same script commenting in/out these 2 lines.

I'm wondering if this is because the model was not trained with these parameters, since I would hope this slowdown doesn't happen for models that are trained with deterministic=True? @rfriedman22 would love to get your input. I'm also going to try training some models with/without the random seed to see if I can figure out what's going on and how we can incorporate a workaround.

Reference Issues/PRs

No issue to reference.

What does this implement/fix? Explain your changes.

In the case of multi-class classification, each prediction task is not entirely separate. A user might want to compute the micro average of the ROC or PR curves, or compute the F1 score (micro or macro averaged). Computing these metrics requires knowledge of all targets and all predictions, at the same time. Currently the implementation does not allow this because it loops over each column of predictions separately, computes performance metrics, and then averages them together at the end.

I added a condition to the compute_score function that allows for this. In the config file, the user specifies the targets as single values (i.e. a column vector). The NN will output K predictions per example, corresponding to the probability that the example belongs to each of K classes.

Performance of the original implementation can be rescued by either (1) one-hot encoding the targets or (2) using implementations of performance metrics that use macro averaging, and specifying them in the config file.

What testing did you do to verify the changes in this PR?

Ran Selene using data where targets is a column vector of integer values 0,...,K-1 and NN architecture outputs K values for each example, corresponding to probabilities of the example belonging to each of the K classes. Wrote custom wrappers of performance metrics and confirmed that the metrics are only called once per epoch, rather than K times.

…ndom suffix to output directory name

Reference Issues/PRs

What does this implement/fix? Explain your changes.

What testing did you do to verify the changes in this PR?

At present, Selene does not save a copy of python file containing the model source code. This seems like an oversight; the output directory contents alone should be enough to use a model trained with Selene.

I've found that in the case of train_model, none of logging handlers ever close (i.e. selene_sdk.train_model.log, selene_sdk.train_model.train.txt, and selene_sdk.train_model.validation.txt). I assume this happens in evaluate_model too but I haven't looked too closely.

This creates problems because I frequently write code that involves creating multiple trainers. Since this happens in the same script, each trainer ends up sharing loggers. If trainer 1 is made before trainer 2, then trainer 2 logs to the trainer 1 log as well as the trainer 2 log.

It would be useful if somewhere within the TrainModel object, the Handlers could be removed from the loggers at some point. To me it seems most useful after calling train_and_validate but it could also be done in an implementation of __del__.

In a related issue, it would be helpful if I could set a flag telling the logger not to write to stderr. Since the logs are getting written elsewhere anyways, it feels redundant to write the logs to stderr, especially since it gets mixed in with my own custom logging.

Reference Issues/PRs

What does this implement/fix? Explain your changes.

One new script (selene_sdk/updating_seqweaver.py) implements a class which constructs and trains new data for seqweaver. selene_sdk/targets/genomic_features.py was modified to add strand specificity when parsing the bed file. I had to change the indexing order in order to maintain the following bed file column order: [chr, start, end, strand, target]. selene_sdk/targets/_genomic_features.pyx was modified to fix the indexing after adding the strand column.

What testing did you do to verify the changes in this PR?

I tested the update_seqweaver.py script, which uses the other scripts as an API, using the original training dataset for Seqweaver (~218 CLIP-seq datasets for 90 RBPs).

Dear Selene/Sei developers, thank you for the colossal work you've undertaken.

I have a question regarding the variant effect prediction functionality of Sei model. I am using the model to calculate the variant effects using gnomAD, there seem to be many mismatches between the REFs of the gnomAD variants and the GRCh38 assembly fasta file. So, my question is, in these cases, does Selene use the REF of the given VCF file or does it use the corresponding NT from the GRCh38 fasta file?