File: test_pbcore_extract_bam_sequence.py

package info (click to toggle)
python-pbcore 1.2.11%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 18,612 kB
  • ctags: 5,336
  • sloc: python: 22,160; xml: 2,667; makefile: 239
file content (63 lines) | stat: -rw-r--r-- 2,227 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

"""
A corner case: reference sequence contains IUPAC ambiguous nucleotide
characters, which we allow.  This test ensures that we can extract the
transcript from a BAM alignment file, specifically when the reverse complement
is needed.
"""

import pbcore.io
import pysam
import tempfile
import unittest
import shutil
import os

sam_str_ = """\
@HD	VN:1.5\tSO:coordinate\tpb:3.0.1
@SQ\tSN:genome\tLN:28\tM5:734d5f3b2859595f4bd87a2fe6b7389b
@RG\tID:1\tPL:PACBIO\tDS:READTYPE=SUBREAD;DeletionQV=dq;DeletionTag=dt;InsertionQV=iq;MergeQV=mq;SubstitutionQV=sq;Ipd=ip;BASECALLERVERSION=2.0.1.0.123678;FRAMERATEHZ=75.000000;BINDINGKIT=foo;SEQUENCINGKIT=bar\tPU:movie1
@PG\tID:bax2bam-0.0.2\tPN:bax2bam\tVN:0.0.2\tDS:bax2bam
movie1/54130/0_10\t16\tgenome\t12\t30\t10=\t*\t0\t-28\tTCTCAGGAGA\t*\tRG:Z:1\tdq:Z:2222'$22'2\tdt:Z:NNNNAGNNGN\tip:B:C,255,2,0,10,22,34,0,2,3,0,16\tiq:Z:(+#1'$#*1&\tmq:Z:&1~51*5&~2\tnp:i:1\tqe:i:10\tqs:i:0\trq:f:0.854\tsn:B:f,-1,-1,-1,-1\tsq:Z:<32<4<<<<3\tzm:i:54130\tAS:i:-3020\tNM:i:134\tcx:i:2"""

fasta_str = """\
>genome
AAAATGASGAGATCARAATGASGAGATC"""
fai_str = """genome  28      8       28      29"""

class TestCase (unittest.TestCase):
    def setUp (self):
        f = open("tst_pbcore.sam", "w")
        f.write(sam_str_)
        f.close()
        # convert to bam using pysam
        sam_in = pysam.AlignmentFile("tst_pbcore.sam", "r")
        bam_out = pysam.AlignmentFile("tst_pbcore.bam", "wb",
            template=sam_in)
        for s in sam_in:
            bam_out.write(s)
        bam_out.close()
        open("tst_pbcore.fa", "w").write(fasta_str)
        open("tst_pbcore.fa.fai", "w").write(fai_str)

    def tearDown (self):
        os.remove("tst_pbcore.sam")
        os.remove("tst_pbcore.bam")
        os.remove("tst_pbcore.fa")
        os.remove("tst_pbcore.fa.fai")

    def test_reverse_complement (self):
        bam_file = pbcore.io.BamReader("tst_pbcore.bam",
            referenceFastaFname="tst_pbcore.fa")
        for aln in bam_file:
            seq = aln.transcript()

if __name__ == "__main__":
    tmp_dir = tempfile.mkdtemp()
    cwd = os.getcwd()
    os.chdir(tmp_dir)
    try:
        unittest.main()
    finally:
        os.chdir(cwd)
        shutil.rmtree(tmp_dir)