build.rna_tools.tools.rna_alignment.rna_alignment.RNAalignment.__init__()

Complexity

Conditions

14

Size

Total Lines	76
Code Lines	59

Duplication

Lines	0
Ratio	0 %

Importance

Changes

0

#!/usr/bin/env python2
    def __init__(self, fn='', fetch=''):
        if fetch:
            import urllib3
            http = urllib3.PoolManager()
            response = http.request('GET', 'http://rfam.xfam.org/family/' +
                                    fetch + '/alignment/stockholm?gzip=1&download=1')
            if not response.status == 200:
                raise RFAMFetchError(
                    "The alignment could not be downloaded. Please check the RFAM id that you requested! (don't put .stk etc in the id)")
            with open(fetch + '.stk.gz', 'wb') as f:
                f.write(response.data)
            with gzip.open(fetch + '.stk.gz', 'rb') as f:
                file_content = f.read()
            with open(fetch + '.stk', 'wb') as f:
                f.write(file_content)
            fn = fetch + '.stk'

        self.fn = fn
        self.lines = open(fn).read().split('\n')
        self.io = AlignIO.read(fn, "stockholm")
        self.ss_cons = self.get_ss_cons()
        self.ss_cons_pk = self.get_ss_cons_pk()
        self.copy_ss_cons_to_all()
        self._ss_cons_std = self.ss_cons
        self.rf = self.get_gc_rf()
        self.shift = self.get_shift_seq_in_align()

        # get all lines not # nor //
        # fix for blocked alignment
        seq_lines = [l for l in self.lines if (
            not l.startswith('#')) and (not l.startswith('//')) and (l)]
        seqs_dict = OrderedDict()
        for seq in seq_lines:
            seq_id, seq_seq = seq.split()
            if seq_id not in seqs_dict:
                seqs_dict[seq_id] = seq_seq
            else:
                seqs_dict[seq_id] += seq_seq
        self.seqs = []
        for seq in seqs_dict:
            self.seqs.append(RNASeq(seq, seqs_dict[seq], ss=self.ss_cons_with_pk_std))

        # this is sick!
        # I create a Cols object to be able to slice alignments
        class Cols:
            def __init__(self, alignment):
                self.alignment = alignment

            def __getitem__(self, i):
                """Return new alignment"""
                if type(i) is list:
                    pass
                else:
                    # collect "new" sequences
                    n_seqs = []
                    for s in self.alignment:
                        new_seq = RNASeq(s.id, s.seq[i], s.ss[i])
                        n_seqs.append(new_seq)

                    # this is not very smart :(
                    # save new seqs to a file
                    # and load it as RNAalignment
                    tf = tempfile.NamedTemporaryFile(delete=False)
                    tf.name += '.stk'
                    with open(tf.name, 'w') as f:
                        f.write('# STOCKHOLM 1.0\n')
                        for s in n_seqs:
                            f.write(' '.join([s.id, s.seq, '\n']))
                        # add ss_cons & //
                        f.write('#=GC SS_cons ' + self.alignment.ss_cons[i] + '\n')
                        if self.alignment.ss_cons_pk:
                            f.write('#=GC SS_cons_pk' + self.alignment.ss_cons_pk[i] + '\n')
                        f.write('#=GC RF ' + self.alignment.rf[i] + '\n')
                        f.write('//\n')
                    return RNAalignment(tf.name)
        self.cols = Cols(self)
        # ^^^^ sick ^^^^^^^^^^^