class Sequence(object):
    # 使用__init__对Sequence进行初始化，Sequence对象具有三个属性
    def __init__(self, name, seq, descr=''):
        self.name = name
        self.seq = seq
        self.descr = descr

    # 给Sequence对象定义一个反向互补方法
    def reverse_complement(self):
        RC = {'A':'T','T':'A','C':'G','G':'C','N':'N',
              'a':'t','t':'a','c':'g','g':'c','n':'n'}
        return Sequence(self.name, ''.join([RC[i] for i in self.seq[::-1]]), self.descr)

    # 定义输出到文件的方法，接受文件句柄输入
    def write_to_fasta_file(self, file_handle):
        def _SeqFormat(seq, chara=80):
            tmp = ''
            for i in range(0,len(seq),chara):
                tmp += (seq[i:(i+chara)]+'\n')
            return tmp
        file_handle.write('>{} {}\n'.format(self.name, self.descr))
        file_handle.write('{}'.format(_SeqFormat(self.seq)))

class FastaReader(object):
    def __init__(self, fastafile):
        self.fasta = fastafile

    # 使用__iter__方法实现对fasta文件的循环解析，其中用到yeild构造了生成器
    def __iter__(self):
        with open(self.fasta) as f:
            seq = None
            for line in f:
                if line.startswith( ">" ):
                    if seq:
                        yield Sequence( name, seq, descr )
                    name, sep, descr = line[1:-1].partition(' ')
                    seq = ""
                else:
                    assert seq is not None, "FASTA file does not start with '>'."
                    seq += line[:-1].encode()
        if seq is not None:
            yield Sequence( name, seq, descr )

# example
for i in FastaReader('ref.fa'):
    # 使用
    i.write_to_fasta_file(file_handle)
    i.seq
    i.name
    i.descr

# 错误的结果
table = open('qc.stat', 'r')
header = [[],] * len(table.readline().strip().split('\t'))
for i in table:
    for x,y in enumerate(i.strip().split('\t')):
        header[x].append(y)
header
# table.close()

# 正确代码
table = open('qc.stat', 'r')
header = [[i,] for i in table.readline().strip().split('\t')]
for i in table:
    for x,y in enumerate(i.strip().split('\t')):
        header[x].append(y)
header
# table.close()

# 错误原因查找

head = '#library\treads_num\treads_len\tbases_num\tGC\tavgQ\tQ20\tQ30'

header = [[], ] * len(head.strip().split('\t'))
print 'error method...'
for i in header:
    print id(i)

header = [[i,] for i in head.strip().split('\t')]
print '\nright method...'
for i in header:
    print id(i)

#初始化定义
seq = {}    #seq = dict()
seq = {'key1':'value1', 'key2':'value2'}

#key不存在时，返回自定义value
seq.get(key, value)

#key不存在时，设置key对应键值为list；key存在时返回seq[key]
seq.setdefault(key,[])

#设置默认字典，读取seq[key]时进行初始化
from collections import defaultdict
seq = defaultdict(list)

# 方法一
list.count('X')    # int

# 方法二
from collections import Counter
Counter(list)['X']    # dict['X'], int

def parse_fasta(fa):
    seq = {}
    with open(fa) as fa:
        for i in fa:
            if i.startswith(">"):
                name = i[1:].partition(' ')[0]
            else:
                seq[name] = seq.get(name, '') + i.strip()
    return seq

def parse_fasta_by_HTSeq(fa):
    import HTSeq
    return {i.name:i.seq for i in HTSeq.FastaReader(fa)}

def format_fasta(seq, len=60):
    tmp = ''
    for i in range((len(seq)//len)+1):
        tmp += '{}\n'.format(seq[(i*60):((i+1)*60)])
    return tmp

#!/usr/bin/env python
# -*- coding:utf-8 -*-

import os
import sys
import gzip
from collections import defaultdict
from collections import Counter

def GFF_parser(file):
    openfile = gzip.open if file.endswith('gz') else open
    with openfile(file) as anno_file:
        for line in anno_file:
            if line.startswith('#'):
                continue
            seqname, source, feature, start, end, score, strand, \
            frame, attributeStr = line.strip().split('\t')
            attribute = [i.strip() for i in attributeStr.split(';') if i.strip() != '']
            attr = {i.partition('=')[0]:i.partition('=')[2] for i in attribute}
            yield seqname, source, feature, int(start), int(end), score, strand, frame, attr

def format_gtf_line(chro, feature, start, end, strand, attr, \
    source=None, score=None, frame=None):
    source = source if source else 'RefSeq'
    score = score if score else '.'
    frame = frame if frame else '.'
    attrstr = ''.join(['{} "{}"; '.format(i[0], i[1]) for i in sorted(attr.items())])
    gtf_line = '{chro}\t{source}\t{feature}\t{start}\t{end}\t{score}\t\
        {strand}\t{frame}\t{attrstr}'.format(**locals())
    return gtf_line

def main(gff, output):
    gtf_out = open(output, 'w')

# RefSeq数据库中GFF文件feature收集，或去除无意义项，或替换
    skip = ['region', 'miRNA', 'cDNA_match', 'repeat_region', 'D_loop', \
            'match', 'origin_of_replication', 'centromere', 'sequence_feature']
    transcript = ['mRNA', 'lnc_RNA', 'primary_transcript', 'tRNA', 'rRNA', \
            'C_gene_segment', 'V_gene_segment', 'SRP_RNA', 'ncRNA', 'transcript']
    biotype = ['lnc_RNA', 'tRNA', 'rRNA', 'SRP_RNA']

    order = []
    gene_info = defaultdict(dict)
    transcript_info = defaultdict(dict)
    check_gene_name = {}    #检查重复gene_name
    check_transcript_name = {}
    for i in GFF_parser(gff):
        chro, source, feature, start, end, score, strand, frame, attr = i

        if feature in skip:
            continue
        elif feature == 'gene':
            gene_id = attr['ID']
            gene_info[gene_id]['name'] = attr['Name']
            gene_info[gene_id]['biotype'] = attr['gene_biotype']
            check_gene_name[gene_id] = attr['Name']
        #使用字典transcript_info储存转录本信息
        elif feature in transcript:
            transcript_id = attr['ID']
            gene_id = attr['Parent'] if 'Parent' in attr else transcript_id
            transcript_name = attr['transcript_id'] if 'transcript_id' \
                in attr else transcript_id
            order.append(transcript_id)
            transcript_info[transcript_id]['gtf_info'] = [chro, start, end, strand]
            transcript_info[transcript_id]['name'] = transcript_name
            transcript_info[transcript_id]['gene'] = gene_id
            # 部分转录本不存在对应gene信息，直接关联转录类型
            if feature in biotype:
                transcript_info[transcript_id]['biotype'] = feature
            check_transcript_name[transcript_id] = transcript_name
        elif feature == 'exon':
            transcript_id = attr['Parent']
            try:
                transcript_info[transcript_id]['exon'][1].append(start)
                transcript_info[transcript_id]['exon'][2].append(end)
            except KeyError:
                transcript_info[transcript_id]['exon'] = [chro, [start,], \
                    [end,], strand]
        elif feature == 'CDS':
            transcript_id = attr['Parent']
            try:
                transcript_info[transcript_id]['CDS'][1].append(start)
                transcript_info[transcript_id]['CDS'][2].append(end)
                transcript_info[transcript_id]['CDS'][4].append(frame)
            except KeyError:
                transcript_info[transcript_id]['CDS'] = [chro, [start,], \
                    [end,], strand, [frame,]]
        else:
            sys.stderr.write('Warning: Annotation Feature "{}" is not collected, '\
                'Unknown error may occur.\n'.format(feature))

# 检查重复gene_name， transcript_names，进行替换
    dup_name = lambda list_in: set([i[0] for i in Counter(list_in).items() if i[1] > 1])
    gene_name = [i[1] for i in check_gene_name.items()]
    dup_gene_name = dup_name(gene_name)

    transcript_name = [i[1] for i in check_transcript_name.items()]
    dup_transcript_name = dup_name(transcript_name)

# 输出GTF文件，order控制输出顺序
    for transcript in order:
        attr = {}
        gene_id = transcript_info[transcript]['gene']
        try:
            gene_name = gene_info[gene_id]['name']
            gene_name = gene_id if gene_name in dup_gene_name else gene_name
        except KeyError:
            gene_name = gene_id
        attr['gene_id'] = gene_name
        transcript_name = transcript_info[transcript]['name']
        transcript_name = transcript if transcript_name in 、
            dup_transcript_name else transcript_name
        attr['transcript_id'] = transcript_name
        try:
            attr['transcript_biotype'] = transcript_info[transcript]['biotype']
        except KeyError:
            attr['transcript_biotype'] = gene_info[gene_id]['biotype']
        try:
            attr['gene_biotype'] =  gene_info[gene_id]['biotype']
        except KeyError:
            attr['gene_biotype'] =  transcript_info[transcript]['biotype']
        chro, start, end, strand = transcript_info[transcript]['gtf_info']
        gtf_out.write('{}\n'.format(format_gtf_line(chro,'transcript',、
            start,end,strand,attr)))
        try:
            exon_feature = transcript_info[transcript]['exon']
        except KeyError:
            exon_feature = transcript_info[transcript]['CDS'][:-1]
        chro, start, end, strand = exon_feature
        for i,j in enumerate(start):
            gtf_out.write('{}\n'.format(format_gtf_line(chro,'exon',、
                start[i],end[i],strand,attr)))
        try:
            cds_feature = transcript_info[transcript]['CDS']
            chro, start, end, strand, frame = cds_feature
            for i,j in enumerate(start):
                gtf_out.write('{}\n'.format(format_gtf_line(chro, 'CDS', start[i], \
                        end[i], strand, attr, frame=frame[i])))
        except KeyError:
            continue
        #部分注释信息不完善，多次使用try... except...输出替代的对应注释信息
    gtf_out.close()