Python SeqIO.QualityIO类代码示例

    def write_record(self, record):
        """Write a single Phd record to the file."""
        assert record.seq, "No sequence present in SeqRecord"
        # This method returns the 'phred_quality' scores or converted
        # 'solexa_quality' scores if present, else raises a value error
        phred_qualities = QualityIO._get_phred_quality(record)
        peak_locations = record.letter_annotations.get("peak_location")
        if len(record.seq) != len(phred_qualities):
            raise ValueError("Number of phd quality scores does not match "
                             "length of sequence")
        if peak_locations:
            if len(record.seq) != len(peak_locations):
                raise ValueError("Number of peak location scores does not "
                                 "match length of sequence")
        if None in phred_qualities:
            raise ValueError("A quality value of None was found")
        if record.description.startswith("%s " % record.id):
            title = record.description
        else:
            title = "%s %s" % (record.id, record.description)
        self.handle.write("BEGIN_SEQUENCE %s\nBEGIN_COMMENT\n"
                          % self.clean(title))
        for annot in [k.lower() for k in Phd.CKEYWORDS]:
            value = None
            if annot == "trim":
                if record.annotations.get("trim"):
                    value = "%s %s %.4f" % record.annotations["trim"]
            elif annot == "trace_peak_area_ratio":
                if record.annotations.get("trace_peak_area_ratio"):
                    value = "%.4f" % record.annotations[
                        "trace_peak_area_ratio"]
            else:
                value = record.annotations.get(annot)
            if value or value == 0:
                self.handle.write("%s: %s\n" % (annot.upper(), value))

        self.handle.write("END_COMMENT\nBEGIN_DNA\n")
        for i, site in enumerate(record.seq):
            if peak_locations:
                self.handle.write("%s %i %i\n" % (
                    site,
                    round(phred_qualities[i]),
                    peak_locations[i])
                )
            else:
                self.handle.write("%s %i\n" % (
                    site,
                    round(phred_qualities[i]))
                )

        self.handle.write("END_DNA\nEND_SEQUENCE\n")

 def test_solexa_to_sanger(self):
     """Mapping check for FASTQ Solexa (-5 to 62) to Sanger (0 to 62)"""
     # The point of this test is the writing code doesn't actually use the
     # solexa_quality_from_phred function directly. For speed it uses a
     # cached dictionary of the mappings.
     seq = "N"*68
     qual = "".join(chr(64+q) for q in range(-5, 63))
     expected_phred = [round(QualityIO.phred_quality_from_solexa(q))
                       for q in range(-5, 63)]
     in_handle = StringIO("@Test\n%s\n+\n%s" % (seq, qual))
     out_handle = StringIO()
     SeqIO.write(SeqIO.parse(in_handle, "fastq-solexa"),
                 out_handle, "fastq-sanger")
     out_handle.seek(0)
     record = SeqIO.read(out_handle, "fastq-sanger")
     self.assertEqual(str(record.seq), seq)
     self.assertEqual(record.letter_annotations["phred_quality"],
                      expected_phred)

 def test_sanger_to_solexa(self):
     """Mapping check for FASTQ Sanger (0 to 93) to Solexa (-5 to 62)"""
     # The point of this test is the writing code doesn't actually use the
     # solexa_quality_from_phred function directly. For speed it uses a
     # cached dictionary of the mappings.
     seq = "N"*94
     qual = "".join(chr(33+q) for q in range(0, 94))
     expected_sol = [min(62, int(round(QualityIO.solexa_quality_from_phred(q))))
                     for q in range(0, 94)]
     in_handle = StringIO("@Test\n%s\n+\n%s" % (seq, qual))
     out_handle = StringIO()
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter("always", BiopythonWarning)
         SeqIO.write(SeqIO.parse(in_handle, "fastq-sanger"),
                     out_handle, "fastq-solexa")
         self.assertTrue(len(w) <= 1, w)
     out_handle.seek(0)
     record = SeqIO.read(out_handle, "fastq-solexa")
     self.assertEqual(str(record.seq), seq)
     self.assertEqual(record.letter_annotations["solexa_quality"],
                      expected_sol)

 def test_sanger_to_solexa(self):
     """Mapping check for FASTQ Sanger (0 to 93) to Solexa (-5 to 62)"""
     #The point of this test is the writing code doesn't actually use the
     #solexa_quality_from_phred function directly. For speed it uses a
     #cached dictionary of the mappings.
     seq = "N"*94
     qual = "".join(chr(33+q) for q in range(0,94))
     expected_sol = [min(62,int(round(QualityIO.solexa_quality_from_phred(q))))
                     for q in range(0,94)]
     in_handle = StringIO("@Test\n%s\n+\n%s" % (seq,qual))
     out_handle = StringIO("")
     #Want to ignore the data loss warning
     #(on Python 2.6 we could check for it!)
     warnings.simplefilter('ignore', BiopythonWarning)
     SeqIO.write(SeqIO.parse(in_handle, "fastq-sanger"),
                 out_handle, "fastq-solexa")
     warnings.filters.pop()
     out_handle.seek(0)
     record = SeqIO.read(out_handle, "fastq-solexa")
     self.assertEqual(str(record.seq), seq)
     self.assertEqual(record.letter_annotations["solexa_quality"],
                      expected_sol)

 def test_phred_quality_from_solexa(self):
     """Mapping check for function phred_quality_from_solexa"""
     self.assertEqual(1, round(QualityIO.phred_quality_from_solexa(-5)))
     self.assertEqual(1, round(QualityIO.phred_quality_from_solexa(-4)))
     self.assertEqual(2, round(QualityIO.phred_quality_from_solexa(-3)))
     self.assertEqual(2, round(QualityIO.phred_quality_from_solexa(-2)))
     self.assertEqual(3, round(QualityIO.phred_quality_from_solexa(-1)))
     self.assertEqual(3, round(QualityIO.phred_quality_from_solexa(0)))
     self.assertEqual(4, round(QualityIO.phred_quality_from_solexa(1)))
     self.assertEqual(4, round(QualityIO.phred_quality_from_solexa(2)))
     self.assertEqual(5, round(QualityIO.phred_quality_from_solexa(3)))
     self.assertEqual(5, round(QualityIO.phred_quality_from_solexa(4)))
     self.assertEqual(6, round(QualityIO.phred_quality_from_solexa(5)))
     self.assertEqual(7, round(QualityIO.phred_quality_from_solexa(6)))
     self.assertEqual(8, round(QualityIO.phred_quality_from_solexa(7)))
     self.assertEqual(9, round(QualityIO.phred_quality_from_solexa(8)))
     self.assertEqual(10, round(QualityIO.phred_quality_from_solexa(9)))
     for i in range(10,100):
         self.assertEqual(i, round(QualityIO.phred_quality_from_solexa(i)))

 def test_solexa_quality_from_phred(self):
     """Mapping check for function solexa_quality_from_phred"""
     self.assertEqual(-5, round(QualityIO.solexa_quality_from_phred(0)))
     self.assertEqual(-5, round(QualityIO.solexa_quality_from_phred(1)))
     self.assertEqual(-2, round(QualityIO.solexa_quality_from_phred(2)))
     self.assertEqual(0, round(QualityIO.solexa_quality_from_phred(3)))
     self.assertEqual(2, round(QualityIO.solexa_quality_from_phred(4)))
     self.assertEqual(3, round(QualityIO.solexa_quality_from_phred(5)))
     self.assertEqual(5, round(QualityIO.solexa_quality_from_phred(6)))
     self.assertEqual(6, round(QualityIO.solexa_quality_from_phred(7)))
     self.assertEqual(7, round(QualityIO.solexa_quality_from_phred(8)))
     self.assertEqual(8, round(QualityIO.solexa_quality_from_phred(9)))
     for i in range(10,100):
         self.assertEqual(i, round(QualityIO.solexa_quality_from_phred(i)))

def compare_record(old, new, truncate=None):
    """Quality aware SeqRecord comparison.

    This will check the mapping between Solexa and PHRED scores.
    It knows to ignore UnknownSeq objects for string matching (i.e. QUAL files).
    """
    if old.id != new.id:
        raise ValueError("'%s' vs '%s' " % (old.id, new.id))
    if old.description != new.description \
    and (old.id+" "+old.description).strip() != new.description:
        raise ValueError("'%s' vs '%s' " % (old.description, new.description))
    if len(old.seq) != len(new.seq):
        raise ValueError("%i vs %i" % (len(old.seq), len(new.seq)))
    if isinstance(old.seq, UnknownSeq) or isinstance(new.seq, UnknownSeq):
        pass
    elif str(old.seq) != str(new.seq):
        if len(old.seq) < 200:
            raise ValueError("'%s' vs '%s'" % (old.seq, new.seq))
        else:
            raise ValueError("'%s...' vs '%s...'" % (old.seq[:100], new.seq[:100]))
    if "phred_quality" in old.letter_annotations \
    and "phred_quality" in new.letter_annotations \
    and old.letter_annotations["phred_quality"] != new.letter_annotations["phred_quality"]:
        if truncate and [min(q,truncate) for q in old.letter_annotations["phred_quality"]] == \
                        [min(q,truncate) for q in new.letter_annotations["phred_quality"]]:
            pass
        else:
            raise ValuerError("Mismatch in phred_quality")
    if "solexa_quality" in old.letter_annotations \
    and "solexa_quality" in new.letter_annotations \
    and old.letter_annotations["solexa_quality"] != new.letter_annotations["solexa_quality"]:
        if truncate and [min(q,truncate) for q in old.letter_annotations["solexa_quality"]] == \
                        [min(q,truncate) for q in new.letter_annotations["solexa_quality"]]:
            pass
        else:
            raise ValueError("Mismatch in phred_quality")
    if "phred_quality" in old.letter_annotations \
    and "solexa_quality" in new.letter_annotations:
        #Mapping from Solexa to PHRED is lossy, but so is PHRED to Solexa.
        #Assume "old" is the original, and "new" has been converted.
        converted = [round(QualityIO.solexa_quality_from_phred(q))
                     for q in old.letter_annotations["phred_quality"]]
        if truncate:
            converted = [min(q,truncate) for q in converted]
        if converted != new.letter_annotations["solexa_quality"]:
            print
            print(old.letter_annotations["phred_quality"])
            print(converted)
            print(new.letter_annotations["solexa_quality"])
            raise ValueError("Mismatch in phred_quality vs solexa_quality")
    if "solexa_quality" in old.letter_annotations \
    and "phred_quality" in new.letter_annotations:
        #Mapping from Solexa to PHRED is lossy, but so is PHRED to Solexa.
        #Assume "old" is the original, and "new" has been converted.
        converted = [round(QualityIO.phred_quality_from_solexa(q))
                     for q in old.letter_annotations["solexa_quality"]]
        if truncate:
            converted = [min(q,truncate) for q in converted]
        if converted != new.letter_annotations["phred_quality"]:
            print(old.letter_annotations["solexa_quality"])
            print(converted)
            print(new.letter_annotations["phred_quality"])
            raise ValueError("Mismatch in solexa_quality vs phred_quality")
    return True