README update

Author: rhiju

README.md

@@ -4,4 +4,80 @@ Utilities to extra mutate-and-map 2D maps from big library sequencing runs
 ## TODO 
 
  - [ ] output to HDF5 file with [Nseq x Nres x Nres] tensor, instead of .txt file
- - [ ] test splitting into chunks.
+ - [ ] provide M2 example
+ - [ ] provide MOHCA example
+ - [ ] convert MATLAB visualization steps into python and provide example notebook.
+
+## Workflow
+
+### Get the counts with `bam_to_m2.py`
+First run `bam_to_m2.py`:
+
+```
+usage: bam_to_m2.py [-h] [-b [BAM ...]] -s SEQUENCES_FASTA [-mq MAP_QUALITY] [--mutdel_cutoff MUTDEL_CUTOFF] [--trim_cutoff TRIM_CUTOFF] [-o OUT_TAG]
+                    [-n CHUNK_SIZE] [--start_idx START_IDX] [--end_idx END_IDX]
+
+Processes BAM file to create 1D and 2D mut/del profiles
+
+options:
+  -h, --help            show this help message and exit
+  -b, --bam [BAM ...]   BAM-formatted text file, Bowtie2
+  -s, --sequences_fasta, --fasta SEQUENCES_FASTA
+                        Fasta file with DNA/RNA sequences used for Bowtie2
+  -mq, --map_quality MAP_QUALITY
+                        minimum Bowtie2 MAPQ to consider read
+  --mutdel_cutoff MUTDEL_CUTOFF
+                        Filter for maximum number of mut/del in read (default 0 means no filter)
+  --trim_cutoff TRIM_CUTOFF
+                        Filter for maximum number of missing terminal residues (trim) in read (default 50; 0 means no filter)
+  -o, --out_tag, --out_file_tag OUT_TAG
+                        Tag for outfile [default: derive from first bamfile]
+  -n, --chunk_size CHUNK_SIZE
+                        split with this number of sequences per chunk
+  --start_idx START_IDX
+                        only do the reference sequences from start_idx onwards [default all]
+  --end_idx END_IDX     only do the reference sequences up to end_idx [default all]
+
+Extremely basic processing of CIGAR string in bowtie2 BAM file to produce alignments and then coding.
+```
+
+The output for each input bam will be 9 files, with suffixes: `counts_del_del.txt.gz counts_del_mut.txt.gz counts_del.txt.gz counts_mutdel_mutdel.txt.gz counts_mut_del.txt.gz counts_mutdel.txt.gz counts_mut_mut.txt.gz counts_mut.txt.gz coverage.txt.gz`
+
+The files are:  
+ - `counts_del_del.txt.gz`, ... `counts_mutdel_mutdel.txt.gz` are 2D files. These hold matrices in plain text of integer counts. The count how many reads have events -- `del` (deletion), `mut` (mutation), or `mutdel` (either deletion or mutation) -- co-occuring at residue i and j.  
+	 - These 2D files take the Nres x Nres matrix for each probed reference sequence and concatenate them. The files have Nres values per row and the number of rows is Nres x Nseq, i.e., Nres rows for each reference sequence.
+ - `counts_del.txt.gz`,...`counts_mut.txt.gz` are 1D files. They hold number of reads in which each event occurs at position i. The files have Nres values per row, and the number of rows is Nseq, one for each sequence.
+ - `coverage.txt.gz` hold the number of reads that contributed counts to each position of each sequence. This file has Nres values per row, and the number of rows is Nseq, one for each sequence.
+
+
+**Tips**:  
+ - Bam file should be sorted!  
+ - If you have a huge file, e.g an Ultima-generated .cram or .bam with 1B reads, you can parallelize by using `--start_idx` and `--end_idx` to select out chunks to run on different processes.  
+ - If you have a bunch of `.bam` files, e.g., in different directories from a pre-split UBR run, you can supply them all with `-b` and the outputs will be in the same directories as the bam file. You can then combine the files with `merge_m2.py` (type `merge_m2.py -h` for help).
+ - Merging a bunch of M2 files can get memory intensive! If not using `--start_idx`/`--end_idx`, you can provide `--chunk_size 100` to get more manageable files, and then after merging the files chunk-wise across multiple directories, concatenate them with `cat [chunk1] [chunk2]`.
+
+ 
+### Visualize the maps
+Once you've got the output of `bam_to_m2.py`, you can fire up MATLAB and run:
+
+```
+% Read in data
+[all_counts_3d, all_counts, all_coverage, tags] = read_m2_output( filedir, data_type );
+
+% derive maps from raw counts. 
+all_cov2d = get_all_cov2d( all_counts_3d, all_counts, all_coverage [, BLANK_OUT5, BLANK_OUT3] );
+
+% show maps with highest coverage first
+coverage = get_coverage( all_coverage );
+colorscale = [-5,5];
+[~,sortidx] = sort(coverage,'descend');
+show_m2_maps( all_cov2d, sortidx, colorscale [, BLANK_OUT5, BLANK_OUT3, coverage, headers, tags] ); 
+```
+
+**Tips**:  
+- For `data_type`, use `mutdel` to capture both mutations and deletions in mutate-and-map experiments correlating SHAPE or DMS signals to intrinsic mutations and deletions or to other SHAPE/DMS events. Use `del` for MOHCA-MaP data, where it seems that only deletions carry signal.  
+- `cov2d` maps are `P(i,j)/P(i)P(j) - 1`, i.e. the 2D signal above the background from 1D events.  
+- An alternative is to use `get_all_m2(...)`. In that case, `colorscale = [0,0.05]` works better in `show_m2_maps`.
+ 
+ 
+ 

bam_to_m2.py

@@ -8,7 +8,7 @@
 import shutil
 
 parser = argparse.ArgumentParser(
-                    prog = 'bam_to_2d.py',
+                    prog = 'bam_to_m2.py',
                     description = 'Processes BAM file to create 1D and 2D mut/del profiles',
                     epilog = 'Extremely basic processing of CIGAR string in bowtie2 BAM file to produce alignments and then coding.')