Merge pull request #226 from golobor/update_docs

update docs

doc/technotes.rst → doc/designnotes.rst

@@ -1,5 +1,5 @@
-Technical notes
-===============
+Design notes
+=============
 
 Designing scientific software and formats requires making a multitude of 
 tantalizing technical decisions and compromises. Often, the reasons behind a 

doc/examples/pairtools_phase_walkthrough.ipynb

@@ -7,8 +7,6 @@
    "source": [
     "# Pairtools phase walkthrough\n",
     "\n",
-    "## Pairtools Phase Walkthrough\n",
-    "\n",
     "Welcome to the pairtools phase walkthrough! This notebook will guide you through the process of resolving contacts between homologous chromosomes using haplotype-resolved Hi-C analysis.\n",
     "\n",
     "## What is haplotype-resolved Hi-C?\n",
@@ -36,50 +34,31 @@
    "source": [
     "Several approaches have been developed to process Hi-C data from haplotype-resolved experiments. In `pairtools`, we implement the approach that was used in Erceg et al. Here is its brief outline:\n",
     "\n",
-    "1. [Create the reference genome](#Create-the-reference-genome): create a \"concatenated\" reference genome that contains sequences of both homologs of each chromosome. \n",
+    "1. Create the haplotype-resolved genome. First, we will create a \"concatenated\" reference genome that contains sequences of both homologs of each chromosome.  \n",
     "\n",
     "    - Incorporate known SNVs (usually in .vcf format) into the reference genome using [bcftools](https://samtools.github.io/bcftools/bcftools.html) to create FASTA files with the sequences of both homologs.\n",
     "    - Add suffixes to the name of each homolog that identify the type (`_hap1` or `_hap2`).\n",
     "\n",
     "2. Map the Hi-C data to the concatenated reference and parse resulting alignment into Hi-C pairs. Compared to the standard Hi-C pipeline, this step would contain a couple of modifications:\n",
-    "    - parse allowing multimappers (mapq 0). \n",
-    "    - make the aligner report two suboptimal alignments (aka the second and the third hit).\n",
+    "    - Make the aligner report two suboptimal alignments (aka the second and the third hit).\n",
+    "    - Parse allowing multimappers (mapq 0). \n",
     "    \n",
     "    Note that, upon mapping to the homolog-resolved genome, Hi-C reads will report the identity of their homologue as the suffix of the chromosome name.\n",
     "    \n",
-    "    See sections:\n",
-    "    \n",
-    "    (i) [Download data](#Download-data)\n",
-    "    \n",
-    "    (ii) [Map data with bwa mem to diploid genome](#Map-data-with-bwa-mem-to-diploid-genome)\n",
-    "    \n",
-    "    (iii) [pairtools parse](#pairtools-parse)\n",
-    "    \n",
+    "3. Phase the resulting pairs based on the reported suboptimal alignments. \n",
     "\n",
-    "3. [pairtools phase](#pairtools-phase): phase the pairs based on the reported suboptimal alignments. \n",
+    "    By checking the scores of two suboptimal alignments, we will distinguish the true multi-mappers from unresolved pairs (i.e. cases when the read aligns to the location with no distinguishing SNV). Phasing will remove the haplotype suffixes from chromosome names and add extra fields to the .pairs file with:\n",
     "\n",
-    "    By checking the scores of two suboptimal alignments, we will distinguish the true multi-mappers from unresolved pairs (i.e. cases when the read aligns to the location with no distinguishing SNV).\n",
-    "    Phasing procedure will remove the haplotype suffixes from chromosome names and add extra fields to the .pairs file with:\n",
-    "    \n",
-    "    '.' (non-resolved)\n",
-    "    \n",
-    "    '0' (first haplotype) or \n",
-    "    \n",
-    "    '1' (second haplotype). \n",
-    "    \n",
+    "    - '.' (non-resolved)\n",
+    "    - '0' (first haplotype)    \n",
+    "    - '1' (second haplotype)\n",
     "    \n",
     "    Phasing schema: \n",
     "    \n",
-    "![image.png](attachment:62e74fba-c1c1-44b5-a3e2-3699c3cac7ce.png)\n",
-    "\n",
+    "    ![image.png](attachment:62e74fba-c1c1-44b5-a3e2-3699c3cac7ce.png)\n",
     "\n",
-    "4. Post-procesing. Sort and dedup Hi-C pairs and calculate stats, similarly to the standard Hi-C pipeline. \n",
     "\n",
-    "    See sections:\n",
-    "    \n",
-    "    (i) [pairtools dedup](#pairtools-dedup)\n",
-    "    \n",
-    "    (ii) [Stats](#Stats)"
+    "4. Post-procesing. Sort and [dedup](#pairtools-dedup) Hi-C pairs and calculate [stats](#Stats), similarly to the standard Hi-C pipeline. "
    ]
   },
   {
@@ -111,7 +90,7 @@
    "id": "5ab026af-fe25-4a70-82ef-52af6fb25371",
    "metadata": {},
    "source": [
-    "## Create the reference genome\n",
+    "## Create the homolog-resolved genome\n",
     "\n",
     "To phase input reads, we need to map the data to the concatenated genome with two haplotypes. \n",
     "Below, we will generate such genome in several steps. You will need a reference genome, and one or two lists of mutations to instroduce to the reference.\n",
@@ -363,15 +342,15 @@
    "id": "dfd7c4cb-31dd-43df-8510-95fd0ff9f78f",
    "metadata": {},
    "source": [
-    "#### Create the index of concatenated haplotypes"
+    "#### Create the bwa index of homolog-resolved genome"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "99d28f6f-b754-4a95-95d5-9e5e51d14571",
    "metadata": {},
    "source": [
-    "Concatenate the genomes and index them together. "
+    "Concatenate the genomes of two homologs and index them together. "
    ]
   },
   {
@@ -588,7 +567,8 @@
    "id": "a8a9fd19",
    "metadata": {},
    "source": [
-    "# Map data with bwa mem to diploid genome\n",
+    "## Map data with bwa mem to diploid genome\n",
+    "\n",
     "In homolog-resolved Hi-C experiments, reads are first aligned against the reference genome and then parsed with pairtools, similar to the standard pairtools-based Hi-C pipeline. However, an additional challenge arises in distinguishing between un-phaseable reads (reads that map equally well to two homologous locations on two homologous chromosomes) and multimappers (reads that map to repeats in the genome).\n",
     "\n",
     "To differentiate between these cases, we can examine the top three candidate alignments for each read:\n",

doc/index.rst

@@ -72,7 +72,8 @@ Contents:
    sorting
    formats
    stats
-   technotes
+   protocols_pipelines
+   designnotes
    cli_tools
 
 .. toctree::

doc/protocols_pipelines.rst

@@ -0,0 +1,118 @@
+Workflows and Parameters
+========================
+
+This page provides guidance on using pairtools for the most common Hi-C protocols and 
+helps users fine-tune the pipeline for different variations of the Hi-C protocol. 
+It covers recommended parameters and best practices for processing Hi-C data using pairtools.
+
+Typical Hi-C Workflow
+----------------------
+
+A typical pairtools workflow for processing standard Hi-C data is outlined below. 
+Please, note that this is a shorter version; you can find a more detailed and reproducible example in chapter :ref:`examples/pairtools_walkthrough`.
+
+1. Align sequences to the reference genome with ``bwa mem``:
+
+    .. code-block:: console
+
+        bwa mem -SP index_file input.R1.fastq input.R2.fastq > input.sam
+
+2. Parse alignments into Hi-C pairs using ``pairtools parse``:
+
+    .. code-block:: console 
+
+        pairtools parse -c /path/to/chrom_sizes -o output.pairs.gz input.sam
+
+3. Sort pairs using ``pairtools sort``:
+
+
+    .. code-block:: console
+
+
+        pairtools sort --nproc 8 -o output.sorted.pairs.gz output.pairs.gz
+
+4. Detect and remove duplicates using ``pairtools dedup`` and generate statistics:
+
+    .. code-block:: console
+
+        pairtools dedup \
+        --output output.nodups.pairs.gz \
+        --output-dups output.dups.pairs.gz \
+        --output-unmapped output.unmapped.pairs.gz 
+        --output-stats output.dedup.stats \
+        output.sorted.pairs.gz
+
+5. Aggregate into a cooler file:
+
+    .. code-block:: console
+
+        cooler cload pairs -c1 2 -p1 3 -c2 4 -p2 5 /path/to/chrom_sizes:1000 output.nodups.pairs.gz output.1000.cool
+
+
+
+
+Recommended pairtools parameters for standard Hi-C protocols
+------------------------------------------------------------
+
+To adapt the standard workflow for common variations of the Hi-C protocol, consider adjusting the following parameters:
+
+1. ``pairtools parse --walks-policy``: 
+
+This parameter determines how pairtools parse handles reads with multiple alignments (walks). We recommend specifying the value explicitly, as the default has changed between versions of ``pairtools parse``.
+
+Our current recommendation is to use ``--walks-policy 5unique``, which is the default setting in the latest version of pairtools. With this option, pairtools parse reports the two 5'-most unique alignments on each side of a paired read as a pair. 
+
+This option increases the number of reported pairs compared to the most conservative ``--walks-policy mask``. However, it's important to note that ``5unique`` can potentially report pairs of non-directly ligated fragments (i.e., two fragments separated by one or more other DNA fragments). Such non-direct (also known as "higher-order" or "nonadjacent") ligations have slightly different statistical properties than direct ligations, as illustrated in several Pore-C papers  [`1 <https://www.biorxiv.org/content/10.1101/833590v1.full>`_ , `2 <https://www.nature.com/articles/s41467-023-36899-x>`_].
+
+An alternative is the ``--walks-policy 3unique`` policy, which reports the two 3'-most unique alignments on each side of 
+    a paired read as a pair, thus decreasing the chance of reporting non-direct ligations. 
+    However, ``3unique`` may not work well in situations where the combined length of a read pair is longer than the length of a DNA fragment (e.g. long read experiments). 
+    In this case, the 3' sides of the two reads will cover the same locations in the DNA molecule, and the 3' alignments may end up identical.
+
+Finally, the experimental ``--walks-policy all`` option reports all alignments of a read pair as separate pairs. This option maximizes the number of reported pairs. The downside is that it breaks the assumption that there is only one pair per read,  which is not compatible with retrieval of .sam records from .pairsam output and may also complicate the interpretation of pair statistics.
+
+2. ``pairtools select "(mapq1>=30) and (mapq2>=30)"``: 
+
+This filtering command selects only pairs with high-quality alignments, 
+   where both reads in a pair have a mapping quality (MAPQ) score of 30 or higher. 
+   Applying this filter helps remove false alignments between partially homologous sequences, which often cause artificial high-frequency interactions in Hi-C maps. 
+   This step is essential for generating maps for high-quality dot calls.
+
+   Note that we recommend storing the most comprehensive, unfiltered list of pairs and applying the filter on the fly prior to contact aggregation:
+
+    .. code-block:: console
+
+        pairtools select "(mapq1>=30) and (mapq2>=30)" output.nodups.pairs.gz | \
+            cooler cload pairs -c1 2 -p1 3 -c2 4 -p2 5 chromsizes.txt:1000 - output.mapq_30.1000.cool
+
+
+Technical tips
+--------------
+
+- **Pipe between commands to save space and I/O throughput**
+
+    Use Unix pipes to connect the output of one command directly to the input of the next command in the pipeline. 
+    This eliminates the need to store intermediate files on disk, saving storage space and reducing I/O overhead.
+    Specifically, mapping, parsing, sorting and deduplication can all be connected into a single pipeline:
+
+    .. code-block:: console
+
+        bwa mem -SP index input.R1.fastq input.R2.fastq | \
+        pairtools parse -c chromsizes.txt | \
+        pairtools sort | \
+            --output output.nodups.pairs.gz \
+            --output-dups output.dups.pairs.gz \
+            --output-unmapped output.unmapped.pairs.gz 
+            --output-stats output.dedup.stats
+
+- **Use recommended compression for efficient storage and processing.** .sam, .pairs and .pairsam files are text-based format that are rather inefficient and slow to process.  
+  Pairtools recognize .bam, .gz and .lz4 file extensions and automatically compress and decompress files on the fly.
+  Compression saves space, and reduces I/O overhead at a relatively minor CPU cost.
+
+- **Parallelize tasks and manage resources effectively for faster execution.**
+  Each pairtool has the CLI flags --nproc-in and --nproc-out to control the number of cores dedicated 
+  to input decompression and output compression. Additionally, `pairtools sort` parallelizes sorting with `--nproc`.ß
+
+Example Workflows
+-----------------
+