import pandas
from pprint import pprint, pformat

subworkflow hits:
    configfile: 'blank.conf'
    snakefile: 'hits.snake'

subworkflow ref_index:
    configfile: 'blank.conf'
    snakefile: 'ref_index.snake'
    workdir: 'hits'

sample_info = pandas.read_csv('sample_info.csv')
samples = sample_info.sample_name

# the taxid 28384 is for "other sequences", a top level taxonomy for
# synthetic, artificial, vector, transposon sequences
taxonids = {
    'cni': '6073',
    'smb': '252141',
    'meta': '33208',
    'euk': '2759',
    'archaea': '2157',
    'fungi': '4751',
    'bacteria': '2',
    'virus': '10239',
    'other_seqs': '28384'
}

taxclasses = list(taxonids.keys())
taxclasses.remove('euk')
taxclasses.remove('meta')
taxclasses.extend(['meta_other','euk_other','org_other','nohit'])

taxclass_fa = expand('taxclass_fa/taxclass_{taxclass}_transcripts.fa',
    taxclass=taxclasses
)

taxclass_counts = expand('all_salmon_quant__{taxclass}.csv',
    taxclass=taxclasses
)
filt_taxclass = [_.replace('all_','filt_') for _ in taxclass_counts]
de_taxclass = expand('{taxclass}_de_site.csv',
    taxclass=taxclasses
)

rule all:
    input:
        'Holobiont_transcripts.fa',
        'Holobiont_cdhit.fa.clstr',
        'Holobiont_transcripts_blast_nr_hits.csv.gz',
        'Holobiont_transcripts_minimap_nt_hits.csv.gz',
        'Holobiont_transcripts_hits.gtf',
        'Holobiont_transcripts_hits_tidmap.csv',
        'Holobiont_transcripts_hits_gene_annotation.csv',
        taxclass_fa,
        filt_taxclass

rule build_taxa_graph:
    input: 'hits/nodes.dmp',
    output: 'subtaxa.json'
    run:
        import csv
        from functools import reduce
        from itertools import combinations
        import json
        import networkx as nx
        import pickle
        import sys

        g = nx.DiGraph()

        with open(input[0],'rt') as f :
            for r in csv.reader(f,delimiter='|') :
                r = [_.strip() for _ in r]
                g.add_edge(r[1],r[0])

        subtaxa = {}
        for k,tid in taxonids.items() :
            subtaxa[k] = nx.algorithms.dag.descendants(g,tid).union(set(tid))

        # logic:
        #   Cnidaria
        #   Symbiodiniaceae
        #   Metazoa_other = Metazoa - Cnidaria
        #   Fungi
        #   Eukaryote_other = Eukaryote - Metazoa - Fungi - Symbiodiniaceae
        #   Bacteria
        #   Archaea
        #   Viruses
        #   Transposon
        #   Other Sequences = (taxid 28384) - Transposon
        #   Org Other = everything not found in above

        subtaxa['meta_other'] = subtaxa['meta'].difference(subtaxa['cni'])
        subtaxa['euk_other'] = subtaxa['euk'].difference(
            subtaxa['meta'].union(subtaxa['fungi']).union(subtaxa['smb'])
        )

        # we aren't interested in the meta and euk high level taxa
        del subtaxa['meta']
        del subtaxa['euk']

        selected_taxa = reduce(lambda a,b: a.union(b), subtaxa.values())
        subtaxa['org_other'] = set(g.nodes).difference(selected_taxa)

        # make sure all of the taxonomic classes are disjoint
        # apparently taxids {1,2,4,5,6,7,9} are in common to many of these
        # remove any taxids from categories that are not disjoint
        # there aren't many
        
        for (ka,a),(kb,b) in combinations(subtaxa.items(),2) :
            print('{} vs {}'.format(ka,kb))
            common = a.intersection(b)
            if len(common) > 0 :
                print('  eliminating: {}'.format(common))
                subtaxa[ka] = subtaxa[ka].difference(common)
                subtaxa[kb] = subtaxa[kb].difference(common)


        taxa = set()
        with open(output[0],'wt') as f :
            json.dump({k:list(v) for k,v in subtaxa.items()},f)

rule cdhit_orig:
    input: 'Holobiont_transcripts.fa'
    output: 'Holobiont_cdhit.fa.clstr'
    params: 'Holobiont_cdhit.fa'
    threads: 28
    log: 'Holobiont_cdhit.log'
    shell:
        # -G use global sequence identity
        # -c sequence identity threshold, default 0.9
        # -aS alignment coverage for the shorter sequence, default 0.0
        # -aL alignment coverage for the longer sequence, default 0.0
        # -g  If set to 1, the program will cluster it into the most similar cluster that meet the threshold
        # -M max available memory (Mbyte), default 400
        # -d length of description in .clstr file, default 20
        # -sc 1 - sort clusters by number of sequences
        # -d 0 - write out full fasta header up to first space in cluster file
        # -r 0 - only compare like strands
        # -n 9 - word length
        '''
        cd-hit-est -i {input} -o {params} -G 0 -c 0.95 -aS 0.9 -aL 0.3 -T {threads} -g 1 -sc 1 -d 0 -M 0 -n 9 -r 0 > {log}
        '''

rule cdhit_clustermap:
    input:
        '{path}.clstr'
    output:
        '{path}.clstr.map.json'
    run:
        import json
        import re

        base_patt = r'(?P<id>[0-9]+)\t(?P<len>[0-9]+)nt, >(?P<name>.*)\.\.\. '
        exem_patt = re.compile(base_patt+r'[*]')
        memb_patt = re.compile(base_patt+r'at (?P<qstart>[0-9]+):(?P<qend>[0-9]+):(?P<sstart>[0-9]+):(?P<send>[0-9]+)/./(?P<pid>.*)')

        #>Cluster 0
        #0	587nt, >Acer_BB1_Locus_17414_Transcript_4of7_Confidence_0.785_Length_587... at 1:587:292:880/+/96.10%
        #29	1040nt, >Acer_BB1_Locus_112322_Transcript_1of1_Confidence_1.000_Length_1040... *

        cluster_map = {}
        def add_to_map(curr_cluster) :
            exemplar = [_ for _ in curr_cluster if _['exemplar']]
            assert len(exemplar) == 1

            exemplar = exemplar[0]

            for seq in curr_cluster :
                seq['parent_tid'] = exemplar['qseqid']
                cluster_map[seq['qseqid']] = seq

        with open(input[0],'rt') as f :
            curr_cluster = []
            for i, line in enumerate(f) :
                if line.startswith('>') :
                    if len(curr_cluster) > 0 :
                        add_to_map(curr_cluster)
                        curr_cluster = []
                else :
                    line = line.strip()
                    if line.endswith('*') :
                        match = exem_patt.match(line)
                        if match is None :
                            print(line)
                            assert False

                        d = {
                            'qseqid': match.group('name'),
                            'qstart': 0,
                            'qend': int(match.group('len')),
                            'sstart': 0,
                            'send': int(match.group('len')),
                            'exemplar': True
                        }
                    else :
                        match = memb_patt.match(line)
                        if match is None :
                            print(line)
                            assert False

                        d = {
                            'qseqid': match.group('name'),
                            'qstart': int(match.group('qstart')),
                            'qend': int(match.group('qend')),
                            'sstart': int(match.group('sstart')),
                            'send': int(match.group('send')),
                            'exemplar': False
                        }

                    curr_cluster.append(d)
            add_to_map(curr_cluster)

        with open(output[0], 'wt') as f :
            json.dump(cluster_map,f,indent=2)

rule blast_to_hits:
    input:
        blast=hits('Holobiont_transcripts_blast_nr.tab.gz'),
        acc=ref_index('acc_db.sqlite3')
    output:
        'Holobiont_transcripts_blast_nr_hits.csv.gz'
    run:
        from annotate_transcripts import hits_fields
        import csv
        import gzip
        import math
        import logging
        import sqlite3
        logging.basicConfig(level=logging.INFO)

        blast_cols = [
            'qseqid', 'sseqid', 'stitle', 'pident', 'length',
            'qstart', 'qend', 'sstart', 'send', 'evalue', 'staxid'
        ]

        con = sqlite3.connect(input.acc)
        i = 0

        with gzip.open(output[0],'wt') as f :
            fout = csv.DictWriter(f,fieldnames=hits_fields)
            with gzip.open(input.blast,'rt') as f :
                for r in csv.DictReader(f,blast_cols,delimiter='\t') :

                    r_acc = {}

                    with con :
                        c = con.execute('SELECT '
                            'GeneID, acc, Symbol '
                            'FROM acc WHERE acc=?',
                            (r['sseqid'],)
                        )
                        res = c.fetchall()
                        if len(res) >= 1 :
                            if len(res) > 1 :
                                logging.warn('{} rows matching protein accession {}, picking first'.format(len(res),r['sseqid']))
                            r_acc = dict(zip(
                                ['GeneID','acc','Symbol'],
                                res[0]
                            ))
                            r_acc = {k:v for k,v in r_acc.items() if v}

                    evalue = float(r['evalue'])
                    if evalue == 0 :
                        score = 255
                    else :
                        score = -math.log10(evalue)
                    fout.writerow({
                        'qseqid': r['qseqid'],
                        'qlen': r['qseqid'].split('_')[-1],
                        'qstart': r['qstart'],
                        'qend': r['qend'],
                        'gene_id': r_acc.get('GeneID',r['sseqid']),
                        'slen': r['length'],
                        'sstart': r['sstart'],
                        'send': r['send'],
                        'gene_desc': r['stitle'],
                        'rna_acc': None,
                        'protein_acc': r_acc.get('acc',r['sseqid']),
                        'gene_symbol': r_acc.get('Symbol',''),
                        'staxids': r['staxid'],
                        'score': score
                    })

                    if i % 1000 == 0 : logging.info(i)
                    i += 1

rule minimap_to_hits:
    input:
        minimap=hits('Holobiont_transcripts_minimap_nt.tab.gz'),
        acc=ref_index('acc_db.sqlite3'),
        gene_info=ref_index('entrez_db.sqlite3')
    output:
        'Holobiont_transcripts_minimap_nt_hits.csv.gz'
    run:
        from annotate_transcripts import hits_fields
        import csv
        import gzip
        import logging
        from pprint import pprint, pformat
        import sqlite3
        logging.basicConfig(level=logging.INFO)

        minimap_cols = [
            'qseqid','qlen','qstart','qend','strand',
            'sseqid','slen','sstart','send','nmatch','nmapped','qual'
        ]

        acc_con = sqlite3.connect(input.acc)

        gene_con = sqlite3.connect(input.gene_info)

        i = 0
        with gzip.open(output[0],'wt') as f :
            fout = csv.DictWriter(f,fieldnames=hits_fields)
            with gzip.open(input.minimap,'rt') as f :
                for r in csv.reader(f,delimiter='\t') :
                    logging.debug(pformat(r))

                    # one of the fields looks like s1:i:[0-9]+
                    score = [_[5:] for _ in r[len(minimap_cols):] if _.startswith('s1:i:')]

                    r = dict(zip(minimap_cols, r[:len(minimap_cols)]))

                    r_acc = {}

                    staxids = None
                    with acc_con :
                        logging.debug('running query')
                        c = acc_con.execute('SELECT '
                            'tax_id, GeneID, acc, Symbol '
                            'FROM acc WHERE acc=?',
                            (r['sseqid'],)
                        )
                        res = c.fetchall()
                        logging.debug('done with query')
                        if len(res) >= 1 :

                            staxids = ';'.join(set([str(_[0]) for _ in res]))

                            if len(res) > 1 :
                                logging.warn('{} rows matching protein accession {}, picking the one with the most information'.format(len(res),r['sseqid']))
                                logging.info('choices: {}'.format(pformat(res)))
                                res = sorted(res,key=lambda x: sum([_ is None for _ in res]))
                                logging.info('picked: {}'.format(pformat(res[0])))

                            r_acc = dict(zip(
                                ['tax_id','GeneID','accession','Symbol'],
                                res[0]
                            ))
                            r_acc = {k:v for k,v in r_acc.items() if v}
                            r_acc['tax_id'] = staxids

                    assert (r['qseqid'].startswith('Acer') or staxids is not None)

                    gene_id = r_acc.get('GeneID',r['sseqid'])
                    gene_desc = None
                    logging.info('gene id for acc={}: {}'.format(r['sseqid'],gene_id))
                    with gene_con :
                        logging.info('getting gene description for gene id: {}'.format(gene_id))
                        # I know this isn't good practice, but using a query template doesn't
                        # return records for some unknown reason, probably due to the column datatype inference
                        # sqlite does
                        c = gene_con.execute('SELECT description '
                            'FROM gene_info WHERE GeneID="{}"'.format(gene_id)
                        )

                        res = c.fetchall()
                        logging.info('done searching for gene description, found {} records'.format(len(res)))

                        if len(res) > 0 :
                            logging.info('found gene info: {}'.format(res))
                            if len(res) > 1 :
                                logging.error('Too many results for gene id {}'.format(gene_id))
                                logging.error(pformat(res))
                                assert False
                            gene_desc = res[0][0]

                    fout.writerow({
                        'qseqid': r['qseqid'],
                        'qlen': r['qlen'],
                        'qstart': r['qstart'],
                        'qend': r['qend'],
                        'gene_id': gene_id,
                        'slen': r['slen'],
                        'sstart': r['sstart'],
                        'send': r['send'],
                        'gene_desc': gene_desc,
                        'rna_acc': r_acc.get('accession',r['sseqid']),
                        'protein_acc': None,
                        'gene_symbol': r_acc.get('Symbol',''),
                        'staxids': staxids or '',
                        'score': 0 if len(score) == 0 else score[0]
                    })

                    if (i % 1000) == 0 : logging.info(i)
                    i += 1


rule map_hits_to_taxclass:
    input:
        fa='Holobiont_transcripts.fa',
        nr='Holobiont_transcripts_blast_nr_hits.csv.gz',
        nt='Holobiont_transcripts_minimap_nt_hits.csv.gz',
        taxa='subtaxa.json',
        clstr='Holobiont_cdhit.fa.clstr.map.json'
    output:
        gtf='Holobiont_transcripts_hits.gtf',
        nohitfa='Holobiont_transcripts_nohit.fa.gz'
    run:
        from Bio import SeqIO
        from collections import defaultdict
        import csv
        import gzip
        import logging
        from pprint import pprint, pformat
        logging.basicConfig(level=logging.INFO)

        # custom method for annotating hits, etc
        import annotate_transcripts as at
        #   annotate_transcripts
        #   clustered_transcript_hits,
        #   hits_to_gtf
        #   nohits_to_fa

        logging.info('loading transcripts')
        transcripts = {}
        for rec in SeqIO.parse(input.fa,'fasta') :
            transcripts[rec.id] = rec

        logging.info('loading taxa tree')
        with open(input.taxa, 'rt') as f :
            taxa = {k:set(v) for k,v in json.load(f).items()}

        logging.info('loading cluster map')
        with open(input.clstr,'rt') as f :
            cluster_map = json.load(f)

        # some hits have more than one taxid that map to different classes
        # prioritize more general classes over more specific ones
        class_priority = (
            'bacteria',
            'fungi',
            'archaea',
            'virus',
            'org_other',
            'other_seqs',
            'euk_other',
            'meta_other',
            'cni',
            'smb'
        )
        class_priority = list(zip(range(len(class_priority)),class_priority))

        logging.info('loading minimap hits')
        # only pull hits from blast nohits
        with gzip.open(input.nt,'rt') as f :
            minimap_hits = [_ for _ in csv.DictReader(f,fieldnames=at.hits_fields)]

        logging.info('mapping minimap hits')
        # map nohit results to nt
        nt_map = at.annotate_transcripts(
            transcripts.keys(),
            minimap_hits,
            taxa,
            class_priority,
            cluster_map,
            'minimap'
        )
        logging.info(pformat(nt_map['stats']))

        logging.info('loading blast hits for minimap nohits')
        blast_hits = []
        with gzip.open(input.nr,'rt') as f :
            for hit in csv.DictReader(f,fieldnames=at.hits_fields) :
                if hit['qseqid'] in nt_map['nohits'] :
                    blast_hits.append(hit)

        logging.info('mapping blast results')
        # map blast results
        nohit_nr_map = at.annotate_transcripts(
            nt_map['nohits'],
            blast_hits,
            taxa,
            class_priority,
            cluster_map,
            'blast'
        )
        logging.info(pformat(nohit_nr_map['stats']))

        for k, hits in nohit_nr_map['hits'].items() :
            qseqids = [_['qseqid'] for _ in hits]
            assert len(qseqids) == len(set(qseqids))

        assert nohit_nr_map['nohits'].isdisjoint(nohit_nr_map['annot'])
        assert nohit_nr_map['nohits'].union(nohit_nr_map['annot']) == nt_map['nohits']
        assert (nohit_nr_map['nohits']
                   .union(nohit_nr_map['annot'])
                   .union(nt_map['annot'])
               ) == set(transcripts)

        logging.info('sanity checks')

        # sanity check make sure we didn't reannotate anything
        assert set(nt_map['annot']).isdisjoint(nohit_nr_map['annot'])

        # make sure nr nohits aren't annotated in nt_map
        assert set(nt_map['annot']).isdisjoint(nohit_nr_map['nohits'])

        logging.info('writing nohit fasta, {} nohit sequences'.format(
            len(nohit_nr_map['nohits']))
        )
        with gzip.open(output.nohitfa,'wt') as f :
            SeqIO.write([transcripts[_] for _ in nohit_nr_map['nohits']],f,'fasta')

        # cluster the remaining nohits
        logging.info('clustering nohits')
        nohit_cluster_hits = at.clustered_transcript_hits(
            nohit_nr_map['nohits'],
            cluster_map
        )

        assert set(nohit_cluster_hits).isdisjoint(set(nt_map['annot']).union(nohit_nr_map['annot']))

        logging.info('merging all hits')
        all_hits = at.merge_hits(
            nt_map['hits'],
            nohit_nr_map['hits'],
            nohit_cluster_hits
        )
        logging.info('total hits: {}'.format(len(all_hits)))
        logging.info('total transcripts in gtf: {}'.format(sum([len(v) for v in all_hits.values()])))
        logging.info('total transcripts in hits: {}'.format(
            sum([len(v) for k,v in all_hits.items() if not k.startswith('Acer')])
        ))
        logging.info('total transcripts with no hits: {}'.format(
            sum([len(v) for k,v in all_hits.items() if k.startswith('Acer')])
        ))

        logging.info('writing hit gtf')
        at.hits_to_gtf(all_hits,output.gtf)

rule sort_fasta:
    input:
        fa='Holobiont_transcripts.fa',
        gtf='Holobiont_transcripts_hits.gtf'
    output:
        taxclass_fa
    run:
        from Bio import SeqIO
        import csv
        import gzip
        import re

        tid_patt = re.compile('transcript_id "([^"]*)"')
        taxclass_patt = re.compile('taxclass "([^"]*)"')
        taxon_patt = re.compile('taxids "([^"]*)"')

        taxclass_map = {}
        with open(input.gtf,'rt') as f :
            for r in csv.reader(f,delimiter='\t') :
                feature = r[2]
                if feature == 'transcript' :
                    tid = tid_patt.search(r[-1]).group(1)
                    taxclass = taxclass_patt.search(r[-1])
                    taxclass = taxclass.group(1) if taxclass is not None else 'nohit'
                    taxa = taxon_patt.search(r[-1])
                    if taxa is not None :
                        taxa = taxa.group(1)
                    else :
                        taxa = ''
                    taxclass_map[tid] = {
                        'taxclass': taxclass,
                        'taxonIDs': taxa
                    }

        with open(input.fa,'rt') as f :

            # open all the file pointers
            fps = {}
            for tc in taxclasses :
                fps[tc] = open('taxclass_fa/taxclass_{}_transcripts.fa'.format(tc),'wt')

            for rec in SeqIO.parse(f,'fasta') :
                taxinfo = taxclass_map[rec.id]
                fps[taxinfo['taxclass']].write(
                    '>{}_{} taxclass={} taxids={}\n{}\n'.format(
                        rec.id,taxinfo['taxclass'],taxinfo['taxclass'],taxinfo['taxonIDs'],rec.seq
                    )
                )

        [_.close() for _ in fps.values()]


rule salmon_index:
    input:
        'Holobiont_transcripts.fa'
    output:
        directory('Holobiont_transcripts__salmon_index')
    threads: 16
    shell:
        'salmon index -p {threads} --index {output[0]} -t {input}'

rule genemap:
    input: 'Holobiont_transcripts_hits.gtf'
    output: 'Holobiont_transcripts_hits_tidmap.csv'
    run:
        import csv
        import re

        patt = re.compile('transcript_id "([^"]*)"')

        with open(output[0], 'wt') as out_f :
            out_f = csv.writer(out_f,delimiter='\t')
            with open(input[0],'rt') as f :
                for r in csv.reader(f,delimiter='\t') :
                    geneid = r[0]
                    feature = r[2]
                    if feature == 'transcript' :
                        tid = patt.search(r[-1]).group(1)
                        out_f.writerow([tid,geneid])

rule gtf_genes_to_csv:
    input: 'Holobiont_transcripts_hits.gtf'
    output: 'Holobiont_transcripts_hits_gene_annotation.csv'
    run:
        import csv

        #WP_096859208.1
        #blast
        #gene
        #394
        #496
        #.
        #.
        #.
        #qlen "309"
        #qstart "1"
        #qend "309"
        #gene_id "WP_096859208.1"
        #slen "103"
        #gene_desc "WP_096859208.1 alpha-D-glucose phosphate-specific phosphoglucomutase [Escherichia coli]"
        #rna_acc ""
        #protein_acc "WP_096859208.1"
        #gene_symbol ""
        #staxids "562"
        #score "48.356547323513816"
        #taxclass "bacteria"

        gene_recs = {}
        fields = (
            'seqid',
            'gene_id',
            'gene_desc',
            'gene_symbol',
            'rna_acc',
            'protein_acc',
            'staxids',
            'score',
            'taxclass',
            'qstart',
            'qend',
            'slen',
        )

        with open(input[0],'rt') as f :
            with open(output[0],'wt') as f_out :
                f_out = csv.DictWriter(f_out,fieldnames=fields,extrasaction='ignore')
                f_out.writeheader()
                for r in csv.reader(f,delimiter='\t') :
                    seqid, source, feature, start, end, score, strand, frame, attrs = r
                    attrs = [_.split(' ',1) for _ in attrs.strip().split(';')]
                    attrs = {_[0]:_[1][1:-1] for _ in attrs if len(_) == 2}

                    if feature == 'gene' :
                        gene_recs[seqid] = attrs
                    else :
                        gene_rec = gene_recs[attrs['gene_id']].copy()
                        gene_rec['seqid'] = attrs['transcript_id']
                        f_out.writerow(gene_rec)

rule salmon_quant:
    input:
        r1='{sample}_R1trim.fq',
        r2='{sample}_R2trim.fq',
        tidmap='Holobiont_transcripts_hits_tidmap.csv',
        index='Holobiont_transcripts__salmon_index'
    output:
        '{sample}__salmon_quant/quant.genes.sf'
    threads: 16
    params:
        out='{sample}__salmon_quant'
    shell:
        'salmon quant -p {threads} -l A -i {input.index} -g {input.tidmap} -1 {input.r1} -2 {input.r2} -o {params.out}'

rule csvgather:
    input:
        lambda w: expand('{sample}__salmon_quant/quant.genes.sf',sample=samples)
    output:
        'all_salmon_quant.csv'
    shell:
        'csvgather -j 0 -f NumReads -d , -t "s:NumReads:{{dir}}:" -t "s:__salmon_quant::" {input} -o {output}'


rule sort_salmon:
    input:
        cnts='all_salmon_quant.csv',
        gtf='Holobiont_transcripts_hits.gtf'
    output:
        taxclass_counts
    shell:
        'python sort_salmon_by_taxclass.py {input.gtf} {input.cnts}'

rule filter:
    input:
        counts='all_salmon_quant__{taxclass}.csv',
        info='sample_info.csv'
    output:
        'filt_salmon_quant__{taxclass}.csv'
    shell:
        'detk-filter "zero(site) == 0" {input.counts} {input.info} -o {output}'

rule de:
    input:
        counts='filt_salmon_quant__{taxclass}.csv',
        info='sample_info.csv'
    output:
        temp('{taxclass}_raw_de_site.csv')
    shell:
        'detk-de deseq2 "counts ~ site[CC]" {input.counts} {input.info} -o {output}'

rule de_annot:
    input:
        de='{taxclass}_raw_de_site.csv',
        gtf='Holobiont_transcripts_all_blast.gtf'
    output:
        '{taxclass}_de_site.csv'
    run:
        import csv
        import pandas
        import re

        de = pandas.read_csv(input.de,index_col=0)

        gene_names = {}
        name_re = re.compile('gene_name "([^"]+)"')
        with open(input.gtf,'rt') as f :
            for r in csv.reader(f,delimiter='\t') :
                if r[2] == 'gene' :
                    gene_names[r[0]] = name_re.search(r[-1]).group(1)

        de_cols = de.columns.tolist()

        de['full_name'] = pandas.Series(gene_names)

        de = de[['full_name']+de_cols]
        de.to_csv(output[0])
