Merge branch 'master' of https://github.com/dereneaton/ipyrad

ipyrad/analysis/digest_genome.py

@@ -6,6 +6,7 @@
 
 import os
 import gzip
+from ..assemble.utils import comp
 
 
 class DigestGenome(object):
@@ -134,67 +135,65 @@ def run(self):
             seq = seq.replace("\n", "").upper()
 
             # digest scaffold into fragments and discard scaff ends
-            bits = ["1{}1".format(i) for i in seq.split(self.re1)]
+            bits = ["1{}1".format(i) for i in seq.split(self.re1)][1:-1]
 
             # digest each fragment into second cut fragment
-            if self.re2:
+            if not self.re2:
+                bits = [(i, seq.index(i), len(i)) for i in bits]
+
+            else:
                 bits1 = bits
                 bits = []
+                pos = 0
                 for fragment in bits1:
                     fbits = fragment.split(self.re2)
-                    for fbit in fbits:
-                        if fbit:
-
-                            # only keep bit if it has both cut sites
-                            if (fbit[0] + fbit[-1]).count("1") == 1:
-                                bits.append(fbit.strip("1"))
-
-            # filter fragments
-            filtered_bits = []
-            for bit in bits:
-                blen = len(bit)
-                if (blen >= self.min_size) and (blen <= self.max_size):
-                    filtered_bits.append(bit)
+                    if len(fbits) > 1:
+                        # remove the 1
+                        fbit = fbits[0][1:]       # 1----2
+                        rbit = fbits[1][:-1]      # 2----1
+
+                        lef = len(fbit)
+                        if (lef > self.min_size) and (lef <= self.max_size):
+                            #pos = seq.index(fbit)
+                            bits.append((fbit, pos, lef))
+
+                        lef = len(rbit)
+                        if (lef > self.min_size) and (lef <= self.max_size):
+                            res = comp(rbit)[::-1]
+                            #pos = seq.index(res)
+                            bits.append((res, pos, lef))
 
             # turn fragments into (paired) reads
-            read1s = []
-            read2s = []
-            for fragment in filtered_bits:
-                r1 = fragment[:self.readlen]
-                r2 = fragment[-self.readlen:]
-                if self.paired:
-                    read1s.append(r1)
-                    read2s.append(r2)
-                else:
-                    read1s.append(r1)
-                    read1s.append(r2)
-
-            # write reads to a file
             fastq_r1s = []
             fastq_r2s = []            
-            for ridx in range(len(read1s)):
-                read1 = read1s[ridx]
+
+            for fragment in bits:
+                fragment, pos, end = fragment
+                r1 = fragment[:self.readlen]
+                r2 = comp(fragment[-self.readlen:])[::-1]
+
+                # write reads to a file
                 for copy in range(self.ncopies):
                     fastq = "@{name}_loc{loc}_rep{copy} 1:N:0:\n{read}\n+\n{qual}"
                     fastq = fastq.format(**{
                         'name': name,
-                        'loc': iloc, 
+                        'loc': iloc,
+                        # 'pos': pos,
+                        # 'end': end,
                         'copy': copy,
-                        'read': read1, 
-                        'qual': "B" * len(read1),
+                        'read': r1, 
+                        'qual': "b" * len(r1),
                     })
                     fastq_r1s.append(fastq)
 
-                if self.paired:
-                    read2 = read2s[ridx]
-                    for copy in range(self.ncopies):
+                    if self.paired:                   
                         fastq = "@{name}_loc{loc}_rep{copy} 2:N:0:\n{read}\n+\n{qual}"
                         fastq = fastq.format(**{
                             'name': name,
                             'loc': iloc, 
                             'copy': copy,
-                            'read': read2,
-                            'qual': "B" * len(read2),
+                            'read': r2,
+                            'qual': "b" * len(r2),
                         })
                         fastq_r2s.append(fastq)
                 iloc += 1

ipyrad/assemble/clustmap.py

@@ -1979,6 +1979,7 @@ def mapping_reads(data, sample, nthreads, altref=False):
         raise IPyradError(error5)
 
 
+
 def check_insert_size(data, sample):
     """
     check mean insert size for this sample and update 
@@ -2068,7 +2069,9 @@ def bedtools_merge(data, sample):
 
         # estimates and updates hackers_only max inner value
         if not data.hackersonly.max_inner_mate_distance:
-            check_insert_size(data, sample)
+            check_insert_size(data, sample)       
+        #cmd2.insert(2, "1000")
+        #cmd2.insert(2, "-d")
         cmd2.insert(2, str(data.hackersonly.max_inner_mate_distance))
         cmd2.insert(2, "-d")
     #else:

ipyrad/core/params.py

@@ -42,7 +42,7 @@ def __init__(self):
             ("demultiplex_on_i7_tags", False),
             ("declone_PCR_duplicates", False),
             ("merge_technical_replicates", True),
-            ("exclude_reference", False),
+            ("exclude_reference", True),
             ("trim_loci_min_sites", 4),
         ])