Python sam.PaddedSAM类代码示例

    def testQueryTooLong(self):
        """
        If the query sequence is longer than the total of the lengths in the
        CIGAR operations, a ValueError must be raised.
        """
        # This test just returns. It used to be possible to reach the
        # "Query ... not fully consumed when parsing CIGAR string."
        # ValueError in sam.py, prior to the fix of
        # https://github.com/acorg/dark-matter/issues/630 but it is not
        # possible to get a CIGAR string that has a different total length
        # from the sequence length through to our code in sam.py because
        # pysam catches the error.  I'm leaving this test here because it
        # documents that the error checked for in sam.py cannot currently
        # be reached and the test may become useful. For now it just returns.
        return
        data = '\n'.join([
            '@SQ SN:ref1 LN:90',
            'query1 0 ref1 1 60 4M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            error = ('^Query TCTAGG not fully consumed when parsing CIGAR '
                     'string\\.')
            assertRaisesRegex(self, ValueError, error, list, ps.queries())

    def testHardClippingInCIGARButQueryNotHardClipped(self):
        """
        As documented in https://github.com/acorg/dark-matter/issues/630 we
        must deal correctly with a case in which the CIGAR string says a
        query is hard-clipped but the query sequence in the SAM file
        actually isn't. This can be due to a prior alignment with a soft clip,
        in which case the full query sequence has to be given before the
        secondary alignment with the hard clip.
        """
        data = '\n'.join([
            '@SQ SN:Chimp-D00220 LN:8',
            '@SQ SN:D-AM494716 LN:8',
            '@SQ SN:D-XXX LN:8',
            '@SQ SN:Chimp-YYY LN:8',
            'query1 0 Chimp-D00220 1 0 3S5M * 0 0 TTTTGGTT 12345678',
            'query1 256 D-AM494716 1 0 3H5M * 0 0 * *',
            'query1 256 D-XXX 1 0 5H3M * 0 0 * *',
            'query1 0 Chimp-YYY 1 0 8M * 0 0 * *',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read1, read2, read3, read4) = list(ps.queries(addAlignment=True))

            self.assertEqual(Read('query1', 'TGGTT---', '45678!!!'), read1)
            self.assertEqual('TTTTGGTT', read1.alignment.query_sequence)

            self.assertEqual(Read('query1/1', 'TGGTT---', '45678!!!'), read2)
            self.assertEqual('TGGTT', read2.alignment.query_sequence)

            self.assertEqual(Read('query1/2', 'GTT-----', '678!!!!!'), read3)
            self.assertEqual('GTT', read3.alignment.query_sequence)

            self.assertEqual(Read('query1/3', 'TTTTGGTT', '12345678'), read4)
            self.assertEqual('TTTTGGTT', read4.alignment.query_sequence)

    def testAllMMatch(self):
        """
        A simple all-'M' match must result in the expected padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 2 60 6M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)

    def testMixedMatch(self):
        """
        A match that is a mix of M, =, and X must result in the expected
        padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)

    def testQuerySoftClipReachesRightEdge(self):
        """
        A match with a soft-clipped region that reaches to the right edge of
        the reference must result in the expected padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 5 60 2M4S * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', '----TCTAGG', '!!!!ZZZZZZ'), read)

    def testQuerySoftClipProtrudesBothSides(self):
        """
        A match with a soft-clipped region that extends to both the left and
        right of the reference must result in the expected padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 4 60 5S5M5S * 0 0 TCTAGGCTGACTAAG ZZZZZZZZZZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', 'TAGGCTGACT', 'ZZZZZZZZZZ'), read)

    def testReferenceDeletion(self):
        """
        An deletion of reference bases must result in the expected padded
        sequence (with Ns inserted for the deleted reference bases).
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 2 60 2M2D4M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', '-TCNNTAGG-', '!ZZ!!ZZZZ!'), read)

    def testMixedMatchSpecificReferenceButNoMatches(self):
        """
        A request for reads aligned against a reference that exists but that
        has no matches must result in an empty list.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:15',
            '@SQ SN:ref2 LN:15',
            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename, referenceIds={'ref2'}))
            self.assertEqual([], list(ps.queries()))

    def testReferencesToStr(self):
        """
        The referencesToStr method must return the expected string.
        """
        data = '\n'.join([
            '@SQ SN:id1 LN:90',
            '@SQ SN:id2 LN:91',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(filename)
            self.assertEqual('id1 (length 90)\nid2 (length 91)',
                             ps.referencesToStr())
            ps.close()

    def testDropSecondary(self):
        """
        Dropping matches flagged as secondary must give the expected result.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
            'query2 256 ref1 2 60 2= * 0 0 TC ZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename, dropSecondary=True))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)

    def testQuerySoftClipProtrudesLeft(self):
        """
        A match with a soft-clipped region that extends to the left of the
        reference must result in the expected padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 2 60 4S2M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', 'AGG-------', 'ZZZ!!!!!!!'), read)

    def testRcNeeded(self):
        """
        A reverse-complimented match (flag = 16) when rcNeeded=True is passed
        must result in the expected (reverse complimented) padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 16 ref1 2 60 6M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(filename)
            (read,) = list(ps.queries(rcNeeded=True))
            self.assertEqual(Read('query1', '-CCTAGA---'), read)
            ps.close()

    def testRcNeeded(self):
        """
        A reverse-complemented match (flag = 16) when rcNeeded=True is passed
        must result in the expected (reverse complemented) padded sequence
        and reversed quality string.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 16 ref1 2 60 6M * 0 0 TCTAGG 123456',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(SAMFilter(filename))
            (read,) = list(ps.queries(rcNeeded=True))
            self.assertEqual(Read('query1', '-CCTAGA---', '!654321!!!'), read)

    def testReferenceSkipAlternateChar(self):
        """
        An skip of reference bases must result in the expected padded
        sequence (with gaps) when a queryInsertionChar is passed.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 2 60 2M2N4M * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(filename)
            (read,) = list(ps.queries(queryInsertionChar='X'))
            self.assertEqual(Read('query1', '-TCXXTAGG-'), read)
            ps.close()

    def testQuerySoftClipProtrudesRight(self):
        """
        A match with a soft-clipped region that extends to the right of
        the reference must result in the expected padded sequence.
        """
        data = '\n'.join([
            '@SQ SN:ref1 LN:10',
            'query1 0 ref1 6 60 2M4S * 0 0 TCTAGG ZZZZZZ',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            ps = PaddedSAM(filename)
            (read,) = list(ps.queries())
            self.assertEqual(Read('query1', '-----TCTAG'), read)
            ps.close()