Merge pull request #13 from treangenlab/main

mergeback

README.md

@@ -1,9 +1,27 @@
 # MethPhaser
+[![Anaconda-Server Badge](https://anaconda.org/bioconda/methphaser/badges/license.svg)](https://anaconda.org/bioconda/methphaser)
+[![Anaconda-Server Badge](https://anaconda.org/bioconda/methphaser/badges/version.svg)](https://anaconda.org/bioconda/methphaser)
+[![Anaconda-Server Badge](https://anaconda.org/bioconda/methphaser/badges/downloads.svg)](https://anaconda.org/bioconda/methphaser)
+## Install:
+ 1. Through mamba or conda: `mamba install -c bioconda methphaser`, `conda install -c bioconda methphaser`
+ 2. Through github download: 
+ 1. Need to put the executable `methphasing` into `~/.bashrc` for `meth_phaser_parallel` to call it. exact command: `export "PATH=/path/to/methphaser/:$PATH"`
+ 2. Other package requirements: check `methphaser.yaml`
 
-Requirement: 
- check methphaser.yaml
-
-
+## Usage
+!!! NOTE: MethPhaser won't work on sex chromosomes because of random X inactivation! Please only consider autosomes. !!!
+
+!!! We recommend using MethPhaser for at least 30X coverage ONT reads since most of the SNV phasing methods (HapCUT2, Whatshap) only work reasonably above that coverage !!!
+### Example 
+Step 1: Run MethPhaser to get block relationship. MethPhaser default ignores the largest unphased region (`-ml -1`), if the input data does not contain any large poorly mapped region like centromere, please use `-ml -2` to not ignore any gap. 
+
+  ./meth_phaser_parallel -b test_data/HLA.R10.haplotagged.bam -r test_data/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr6.fna -g test_data/LSK.filtered.gtf -vc test_data/HLA.R10.phased.vcf.gz -o test_data/work 
+
+Step 2: Run post processing script to get modified reads and vcf file
+
+  ./meth_phaser_post_processing -ib test_data/HLA.R10.haplotagged.bam -if test_data/work/ -ov test_data/output.vcf -ob test_data/output -vc test_data/HLA.R10.phased.vcf.gz 
+
+```
  usage: meth_phaser_parallel [-h] -b -r -g -vc [-vt] [-t] [-ml] [-c] [-a] [-o] [-k]
 
  methphaser: phase reads based on methlytion informaiton
@@ -25,26 +43,13 @@ Requirement:
  -r , --reference reference genome
  -g , --gtf gtf file from whatshap visualization
  -vc , --vcf_called called vcf file from HapCUT2
-  
-Recomanded phasing flow: 
+```
 
- ![flowchart drawio](https://github.com/treangenlab/methphaser/assets/13065758/c74e5a1d-1c24-49c0-abe3-5b125f43f7eb)
 
-
-## Example 
-Step 1: Run MethPhaser to get block relationship 
-
-  ./meth_phaser_parallel -b test_data/HLA.R10.haplotagged.bam -r test_data/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr6.fna -g test_data/LSK.filtered.gtf -vc test_data/HLA.R10.phased.vcf.gz -o test_data/work -k -1 -ml -2
-
-Step 2: Run post processing script to get modified reads and vcf file
-
-  ./meth_phaser_post_processing -ib test_data/HLA.R10.haplotagged.bam -if test_data/work/ -ov test_data/output.vcf -ob test_data/output -vc test_data/HLA.R10.phased.vcf.gz 
-  
-
-Known Issue:
-1. Does not phase reads before the first phased block and after the last phased block
-2. Multiprocessing threads should be smaller than the total block number on the chromosome. e.g., MethPhaser cannot use 3 threads to phase 3 blocks, because the number of the gaps among 3 blocks is 2. 
-3. MethPhaser sometimes ignore some very short reads (about 0.1% reads). Probably because of pysam. 
+Note: The BAM file should only contain primary alignment (no secondary alignment and supplementary alignment), otherwise it will trigger pysam MM tag reading error. Suggested command:
+`samtools view -bF 2304 -o output.bam input.bam`
 
 
+Recomanded phasing flow: 
 
+ ![flowchart drawio](https://github.com/treangenlab/methphaser/assets/13065758/c74e5a1d-1c24-49c0-abe3-5b125f43f7eb)