Merge pull request #14 from maxibor/dev

CLI update and calmd troubleshooting

docs/source/index.rst

@@ -19,6 +19,7 @@ __ homepage_
    API
    CLI
    output
+   troubleshooting
 
 
 

docs/source/troubleshooting.md

@@ -0,0 +1,11 @@
+# Troubleshooting
+
+## My alignment files don't have a MD tag
+
+You can use [samtools calmd](http://www.htslib.org/doc/samtools-calmd.html) to set the MD tag
+
+Example:
+
+```bash
+samtools calmd -b alignment.bam reference.fasta > aln.bam
+```

pydamage/cli.py

@@ -14,19 +14,16 @@
     default=35,
     type=int,
     show_default=True,
-    help="Window length for damage modeling",
+    help="Window length (in bp) for damage modeling",
 )
 @click.option(
     "-p",
     "--process",
     default=2,
     type=int,
     show_default=True,
-    help="Number of processes",
+    help="Number of processes for parallel computing",
 )
-@click.option("-s", "--show_al", is_flag=True, help="Show alignments representations")
-@click.option("-pl", "--plot", is_flag=True, help="Make the damage plots")
-@click.option("--verbose", is_flag=True, help="Verbose mode")
 @click.option(
     "-o",
     "--outdir",
@@ -35,16 +32,28 @@
     show_default=True,
     help="Output directory",
 )
-@click.option("--force", is_flag=True, help="Force overwriting of results directory")
-@click.option("--group", is_flag=True, help="Group references together for analyis")
+@click.option("-s", "--show_al",
+              is_flag=True,
+              help="Display alignments representations")
+@click.option("-pl", "--plot",
+              is_flag=True,
+              help="Write damage fitting plots to disk")
+@click.option("-vv", "--verbose", is_flag=True, help="Verbose mode")
+@click.option("-f", "--force",
+              is_flag=True,
+              help="Force overwriting of results directory")
+@click.option("-g", "--group",
+              is_flag=True,
+              help="Use entire BAM file as single reference for analyis "
+              "(ignore reference headers)")
 def cli(no_args_is_help=True, **kwargs):
     """\b
     PyDamage: Damage parameter estimation for ancient DNA
     Author: Maxime Borry
     Contact: <borry[at]shh.mpg.de>
     Homepage & Documentation: github.com/maxibor/pydamage
 
-    BAM: path to BAM/SAM/CRAM alignment file
+    BAM: path to BAM/SAM/CRAM alignment file. MD tags need to be set.
     """
     analyze(**kwargs)
 

pydamage/damage.py

@@ -16,7 +16,7 @@ def __init__(self, reference, al_handle):
             al_handle(pysam.AlignmentFile)
         """
         self.alignments = al_handle.fetch(reference)
-
+        self.reference = reference
         # self.alignments = al_file
 
     # def __repr__(self):
@@ -45,14 +45,13 @@ def get_damage(self, wlen, show_al):
         all_bases = []
         for al in self.alignments:
             if al.is_reverse is False and al.is_unmapped is False:
-                cigar = al.cigartuples
-                ref = al.get_reference_sequence()
-                quer = al.query_sequence
 
                 all_damage = damage_al(
-                    reference=ref,
-                    query=quer,
-                    cigartuple=cigar,
+                    reference=al.get_reference_sequence(),
+                    read_name=al.query_name,
+                    ref_name=self.reference,
+                    query=al.query_sequence,
+                    cigartuple=al.cigartuples,
                     wlen=wlen,
                     show_al=show_al,
                 )

pydamage/exceptions.py

@@ -0,0 +1,15 @@
+class Error(Exception):
+    """Base class for exceptions in this module."""
+    pass
+
+
+class AlignmentFileError(Error):
+    """Exception raised for errors in validating the dataset agains the standards.
+
+    Attributes:
+        message -- explanation of the error
+    """
+
+    def __init__(self, message):
+        self.message = message
+

pydamage/main.py

@@ -5,6 +5,7 @@
 from functools import partial
 from pydamage import damage
 from pydamage.plot import damageplot
+from pydamage.exceptions import AlignmentFileError
 from pydamage.accuracy_model import prepare_data, load_model, fit_model
 import sys
 from tqdm import tqdm
@@ -122,6 +123,8 @@ def analyze_multi(
     filt_res = [i for i in res if i]
 
     print(f"{len(filt_res)} contigs were successfully analyzed by Pydamage")
+    if len(filt_res) == 0 :
+        raise AlignmentFileError("Check your alignment file")
 
     if plot and len(filt_res) > 0:
         print("\nGenerating Pydamage plots")
@@ -224,6 +227,10 @@ def analyze_group(
         reflen += i[8]
     cov = cov / nb_ref
 
+    if nb_reads_aligned == 0:
+        raise AlignmentFileError("No Alignments were found\nCheck your alignment file")
+
+
     damage_dict = damage.test_damage_group(
         ct_data,
         ga_data,

pydamage/parse_damage.py

@@ -1,11 +1,13 @@
 #!/usr/bin/env python3
 
 
-def damage_al(reference, query, cigartuple, wlen, show_al):
+def damage_al(reference, read_name, ref_name, query, cigartuple, wlen, show_al):
     """Compute CtoT mutations for a single alignment
 
     Args:
-        reference (string): reference sequence
+        reference (str): reference sequence
+        read_name(str): name of read
+        ref_name(str): name of reference
         query (string): query sequence
         cigartuple (tuple): cigar tuple (pysam)
         wlen (int): window length
@@ -71,13 +73,13 @@ def damage_al(reference, query, cigartuple, wlen, show_al):
                 base_trans_counts['GA'].append(i)
 
     if show_al:
-        res += "R " + r_string + "\n  "
+        res += "R " + r_string + "\t" + ref_name + "\n  "
         for i in range(0, min(len(r_string), len(q_string))):
             if r_string[i] == q_string[i]:
                 res += ("|")
             else:
                 res += (" ")
-        res += "\nQ " + q_string
+        res += "\nQ " + q_string + "\t" + read_name + "\n "
         print(res)
 
     return(base_trans_counts)