Add demultiplexStat for 10X data

workflows/illumina_10X/components/demultiplexstats_10X.py

+#
+# Copyright (C) 2012 INRA
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+from weaver.function import PythonFunction
+import logging
+
+from ng6.analysis import Analysis
+from operator import __getitem__
+
+def write_indices_stats(reads_file, stat_file, state , expected_index_seq = None):
+    import re
+    import jflow.seqio as seqio
+
+    indices_stat = {}
+    total_number = 0
+    total_passing_filter = 0
+    # Count indices
+    reader = seqio.SequenceReader(reads_file)
+    for id, desc, seq, qualities in reader:
+        match = re.search("^\d:[Y|N]:\d+:([ATGCN]+)\+*([ATGCN]*)", desc)
+        print(match)
+        if match:
+            index_seq = match.group(1) + match.group(2)
+            if index_seq not in indices_stat :
+                indices_stat[index_seq] = {}
+                indices_stat[index_seq]["number"] = 0
+                indices_stat[index_seq]["passing_filter"] = 0
+            indices_stat[index_seq]["number"] += 1
+            total_number += 1
+            if re.match("^\d:N:\d+:[ATGCN]+", desc) :
+                indices_stat[index_seq]["passing_filter"] += 1
+                total_passing_filter += 1
+        else:
+            raise ValueError("The description '" + desc + "' of the sequence " + id + " is in an invalid format.")
+
+    if state == 'determined' and len(list(indices_stat.keys())) > 1  and expected_index_seq == None :
+        raise Exception("Please provide the expected index sequence for the file '%s'"%reads_file)
+
+    # expected index seq for mismatch (only for determined)
+    if state == 'determined' and expected_index_seq :
+        indices_stat.clear()
+        indices_stat[expected_index_seq] = {'number' : total_number, "passing_filter" : total_passing_filter }
+
+    # Write indices count
+    fh_stat_file = open(stat_file, 'w')
+    for index_seq in indices_stat.keys():
+       fh_stat_file.write( index_seq + ";" + str(indices_stat[index_seq]["number"]) + ";" + str(indices_stat[index_seq]["passing_filter"]) + "\n" )
+    fh_stat_file.close()
+
+class Demultiplex10XStats (Analysis):
+
+    def define_parameters(self, determined_R1_files, undetermined_R1_files, expected_indexes = None, index_count_threshold=0.8, undetermined_threshold=0.8):
+        self.add_input_file_list( "determined_R1_files", "determined_R1_files", default=determined_R1_files, required=True, file_format = 'fastq')
+        self.add_input_file_list( "undetermined_R1_files", "undetermined_R1_files", default=undetermined_R1_files, required=True, file_format = 'fastq')
+        #Unknown indices with a number of fragments < index_count_threshold*number_of_fragments_in_sample_with_littlest_population are merged in "All others".
+        self.add_parameter("index_count_threshold", "index_count_threshold", default=index_count_threshold, type='float')
+        self.add_parameter("undetermined_threshold", "The percentage of undetermined indexes", default=undetermined_threshold, type='float')
+        self.add_parameter_list('expected_indexes', 'list of expected index for each determined input file, this is used to ensure it allows some mismatches', default = expected_indexes)
+        self.add_output_file_list( "determined_idx_count_files", "determined_idx_count_files", pattern='{basename_woext}.stdout', items=self.determined_R1_files)
+        self.add_output_file_list( "undetermined_idx_count_files", "undetermined_idx_count_files", pattern='{basename_woext}.stdout', items=self.undetermined_R1_files)
+
+    def define_analysis(self):
+        self.name = "Demultiplex10XStats"
+        self.description = "Statistics about 10X demultiplexing"
+        self.software = "-"
+        self.options = str(self.index_count_threshold)
+
+    def post_process(self):
+        # Process samples
+        min_determined = -1
+        indices_stat = self._merged_indices_stats(self.determined_idx_count_files)
+        for index_seq in indices_stat.keys():
+            if min_determined == -1 or min_determined > indices_stat[index_seq]["number"] :
+                min_determined = indices_stat[index_seq]["number"]
+            self._add_result_element(index_seq, "number", str(indices_stat[index_seq]["number"]), "determined")
+            self._add_result_element(index_seq, "passing_filter", str(indices_stat[index_seq]["passing_filter"]), "determined")
+
+        # Process unknown indices
+        other = {"number":0, "passing_filter":0}
+        indices_stat = self._merged_indices_stats(self.undetermined_idx_count_files)
+
+        # check undetermined
+        overmin = 0
+        for data in indices_stat.values():
+            if data["number"] >= min_determined :
+                overmin += 1
+
+        # determine the maximum number of undetermined index (with too much sequences) that have to be saved like new indexs
+        max_nbindexsaved = float(len(list(indices_stat.values()))) # maximum number of undetermined index saved as new indexs
+        if max_nbindexsaved > 100 :
+            max_nbindexsaved = 100
+
+        # Sort undetermined index on number of sequences
+        indices_stat_sorted = sorted(indices_stat, key=lambda x: indices_stat[x]['number'], reverse=True)
+        nbindexsaved = 0
+        for index_seq in indices_stat_sorted:
+            if indices_stat[index_seq]["number"] >= self.index_count_threshold*min_determined and nbindexsaved <= max_nbindexsaved :
+                self._add_result_element(index_seq, "number", str(indices_stat[index_seq]["number"]), "undetermined")
+                self._add_result_element(index_seq, "passing_filter", str(indices_stat[index_seq]["passing_filter"]), "undetermined")
+                nbindexsaved = nbindexsaved + 1
+            else:
+                other["number"] += indices_stat[index_seq]["number"]
+                other["passing_filter"] += indices_stat[index_seq]["passing_filter"]
+        self._add_result_element("All others", "number", str(other["number"]), "undetermined")
+        self._add_result_element("All others", "passing_filter", str(other["passing_filter"]), "undetermined")
+
+
+    def get_version(self):
+        return "-"
+
+    def _merged_indices_stats(self, files):
+        indices_stat = {}
+        for current_file in files:
+            fh_current_file = open(current_file, 'r')
+            for line in fh_current_file:
+                line = line.rstrip()
+                index, number, passing_filter = line.split(';')
+                if index not in indices_stat :
+                    indices_stat[index] = {}
+                    indices_stat[index]["number"] = 0
+                    indices_stat[index]["passing_filter"] = 0
+                indices_stat[index]["number"] += int(number)
+                indices_stat[index]["passing_filter"] += int(passing_filter)
+            fh_current_file.close()
+        return indices_stat
+
+    def process(self):
+        demultiplex_stats = PythonFunction(write_indices_stats, cmd_format="{EXE} {IN} {OUT} {ARG}")
+
+        # determined
+        for idx, infile in enumerate(self.determined_R1_files) :
+            logging.getLogger("ng6").debug("demultiplex10X Stat :")
+            logging.getLogger("ng6").debug(self.expected_indexes)
+            demultiplex_stats( inputs = infile, outputs = self.determined_idx_count_files[idx], arguments= [ 'determined', self.expected_indexes[idx]])
+
+
+        # undetermined
+        for idx, infile in enumerate(self.undetermined_R1_files) :
+            demultiplex_stats( inputs = infile, outputs = self.undetermined_idx_count_files[idx], arguments = 'undetermined')
+
+
+