|
| 1 | +# Find regions of accessible chromatin |
| 2 | + |
| 3 | +Nanopore sequencing can detect multiple base modifications simultaneously and we can leverage this capability by introducing exogenous base modifications at specific functional regions. |
| 4 | +One such method uses a 6mA methyltransferase such as EcoGII or Hia5 to label accessible regions of chromatinized DNA, usually by treatment of cell nuclei with the enzyme. |
| 5 | + |
| 6 | + |
| 7 | + |
| 8 | +## Predict regions of open chromatin |
| 9 | +Modkit comes with a machine learning model that has been trained to identify regions of open chromatin based on 6mA signal. |
| 10 | +You can invoke this model with the following command: |
| 11 | + |
| 12 | +```bash |
| 13 | +$ modkit open-chromatin predict ${mod_bam} \ |
| 14 | + --model ${model} \ |
| 15 | + --log modkit_predict.log \ |
| 16 | + -o ./accessible_regions.bedgraph \ |
| 17 | + --device 0 |
| 18 | +``` |
| 19 | + |
| 20 | +Where `${model}` is the path to the directory with the model, for example `dist_modkit_v0.5.0_38fda16/models/r1041_e82_400bps_hac_v5.2.0@v0.1.0`. |
| 21 | +The output of the command is a bedGraph file with the following schema: |
| 22 | + |
| 23 | +| column | name | description | type | |
| 24 | +|--------|------|-------------|------| |
| 25 | +| 1 | chrom | contig or scaffold | string | |
| 26 | +| 2 | start | start of the region | int | |
| 27 | +| 3 | end | end of the region | int | |
| 28 | +| 4 | probability | accessibility probability, (0.0, 1) | float | |
| 29 | + |
| 30 | +## Configuring `--step-size` |
| 31 | + |
| 32 | +The open chromatin model in Modkit works by making predictions on 100 base pair windows of the genome. |
| 33 | +The output of the model is a prediction _probability_ between 0.0 and 1.0, non-inclusive. |
| 34 | +To make predictions on the whole genome (or regions when the `--include-bed` option is provided) Modkit applies the model to overlapping windows. |
| 35 | +The `--step-size` determines how much to advance before making another prediction: |
| 36 | + |
| 37 | +```text |
| 38 | +Genome <--//------------------//-> |
| 39 | +Window 1 |---w----| |
| 40 | +Window 2 |_s_|---w----| |
| 41 | +``` |
| 42 | +In the above `w` is the window size (100bp with the current model) and `s` is the step size which can be configured via the `--step-size` parameter. |
| 43 | +The smaller the step size, the finer-graned resolution the output at the cost of more computation. |
| 44 | +This can be seen in the following example browser image: |
| 45 | + |
| 46 | + |
| 47 | + |
| 48 | +Step size of 25 base pairs is the default. |
| 49 | + |
| 50 | +## Using `--threshold` to get regions of high confidence |
| 51 | + |
| 52 | +As seen above, the output from `ope-chromatin predict` is a stream of probabilities over the genome or desired region of the genome. |
| 53 | +A lot of these predictions are going to be very close to zero. |
| 54 | +You can remove low-probability regions with the `--threshold` option so that only intervals with a probability greater than or equal to this value are reported. |
| 55 | +Using this option in combination with `bedtools merge` can transform the bedGraph file into a BED of predicted open chromatin regions. |
| 56 | + |
| 57 | +```bash |
| 58 | +$ modkit open-chromatin predict ${bam} \ |
| 59 | + --model ${model} \ |
| 60 | + --threshold 0.8 \ |
| 61 | + --include-bed promoters_slop2000.bed \ |
| 62 | + -o stdout \ |
| 63 | + | bedtools merge -i - > accessible_regions.bed |
| 64 | +``` |
| 65 | + |
| 66 | +## Running with a GPU |
| 67 | + |
| 68 | +Modkit comes distributed with the capability to run open chromatin prediction on a normal CPU, NVIDIA GPUs, and can be built to run on Apple GPUs. |
| 69 | +The normal distribution only has the capability to run on CPU hardware. |
| 70 | +At the current stage of development this configuration is probably too slow for most practical purposes outside of small region checks. |
| 71 | +The `candle` distribution can be downloaded and used with NVIDIA GPUs directly, however we cannot guarantee that it will work on _every_ GPU setup. |
| 72 | +For the best performance on GPU _or_ CPU, the `tch` (pytorch) distribution is recommended, however it requires that you download libtorch from the internet. |
| 73 | +Instructions can be found in the `BUILD_NOTES_*.txt` packaged with the software. |
| 74 | + |
| 75 | +## Troubleshooting: Checking your input data |
| 76 | + |
| 77 | +To predict regions of open chromatin, the input modBAM should have 6mA base modification calls. |
| 78 | +You can quickly check that you have valid 6mA calls using the [`check-tags` command](./intro_modbam_check_tags.md): |
| 79 | + |
| 80 | +```bash |
| 81 | +$ modkit modbam check-tags ${bam} --head 1000 |
| 82 | +``` |
| 83 | + |
| 84 | +This command will return non-zero when you have invalid MM/ML [SAMTags](https://samtools.github.io/hts-specs/SAMtags.pdf). |
| 85 | +You should expect to see an output similar to this: |
| 86 | + |
| 87 | +```text |
| 88 | +> checking tags on first 1000 reads |
| 89 | +> no errors |
| 90 | +> num PASS records: 1000 (100.00%) |
| 91 | +> num records: 1000 |
| 92 | +> valid record tag headers: |
| 93 | ++------------+-------+ |
| 94 | +| tag_header | count | |
| 95 | ++------------+-------+ |
| 96 | +| A+a. | 1000 | |
| 97 | ++------------+-------+ |
| 98 | +
|
| 99 | +
|
| 100 | +> modified bases: |
| 101 | ++--------+--------------+----------+------+ |
| 102 | +| strand | primary_base | mod_code | mode | |
| 103 | ++--------+--------------+----------+------+ |
| 104 | +| + | A | a | . | |
| 105 | ++--------+--------------+----------+------+ |
| 106 | +``` |
| 107 | + |
| 108 | + |
| 109 | + |
0 commit comments