def mode_2(exp_matrix):
    
    #remember value of bedgraph, ugly way
    value = {}
    for regions in exp_matrix.get_regionsets():
        for region in regions:
            value[(region.chrom, region.initial, region.final)] = region.data
    
    for region in exp_matrix.get_regionsets():
        f = open("region_" + str(region.name) + ".data", 'w')
        
        region_set = GenomicRegionSet("")
        _, _, mappedGenes, _, gene_peaks_mapping = region_set.filter_by_gene_association(region.fileName, None, gene_file, genome_file, threshDist=2000)
        
        for k in gene_peaks_mapping.keys():
            chr, raw_positions = k.split(':')
            start, end = map(lambda x: int(x), raw_positions.split('-'))
            
            #if peak is not assigned, an empty string occurs
            if "" in gene_peaks_mapping[k]:
                gene_peaks_mapping[k].remove("")
            
            list = 'NA' if not gene_peaks_mapping[k] else ','.join(gene_peaks_mapping[k])
            
            print(chr, start, end, value[(chr, start, end)], list, sep='\t', file = f)
        
        f.close()

    def get_biotypes(self, gene_set=None):
        """Get the region sets of different Biotypes.

        *Keyword arguments:*

        *Return:*

            - result_grs -- A list of GenomicRegionSets containing the regions for each Biotype.
        """

        # Fetching gene names
        mapped_gene_list = None
        unmapped_gene_list = None

        # Fetching exons
        query_dictionary = {self.GeneField.FEATURE_TYPE: "exon"}
        query_annset = self.get(query_dictionary)

        # Creating GenomicRegionSet
        result_grs = GenomicRegionSet("exon")
        for e in query_annset.gene_list:
            gr = e[self.GeneField.GENOMIC_REGION]
            gr.name = e[self.GeneField.TRANSCRIPT_ID]
            result_grs.add(gr)

        if gene_set:
            return result_grs, unmapped_gene_list
        else:
            return result_grs

def mode_3(exp_matrix):
    #remember value of bedgraph, ugly way
    score = {}
    for regions in exp_matrix.get_regionsets():
        for region in regions:
            score[(region.chrom + ':' + str(region.initial) + '-' + str(region.final))] = region.data
    
    
    for region in exp_matrix.get_regionsets():
        f = open("region_" + str(region.name) + ".data", 'w')
        
        region_set = GenomicRegionSet("")
        _, _, mappedGenes, _, gene_peaks_mapping = region_set.filter_by_gene_association(region.fileName, None, gene_file, genome_file, threshDist=2000)
        
        avg_score = {} #score per peak
        genes = {}
        
        for peak, gene_list in gene_peaks_mapping.items():
            for gen in gene_list: #reverse mapping peak -> gene to gene -> peak
                if not gen:
                    continue
                genes[gen] = genes.get(gen, set())
                genes[gen].add(peak)
                avg_score[gen] = avg_score.get(gen, [])
                avg_score[gen].append(score[peak]) #join all scores of peaks assigned to a gen
        
        for gen in genes.keys():
            avg = sum(map(lambda x: float(x), avg_score[gen]))/ float(len(avg_score[gen]))
            print(gen, avg, ", ".join(str(t) for t in genes[gen]), sep='\t', file = f)
               
        f.close()

    def get_genes(self, gene_set = None):
        """
        Gets regions of genes.
        It returns a GenomicRegionSet with such genes. The id of each gene will be put
        in the NAME field of each GenomicRegion.

        Keyword arguments:
        gene_set -- A set of genes to narrow the search.

        Return:
        result_grs -- A GenomicRegionSet containing the genes.
        unmapped_gene_list -- A list of genes that could not be mapped to an ENSEMBL ID.
        """

        # Fetching gene names
        mapped_gene_list = None
        unmapped_gene_list = None
        if(gene_set): mapped_gene_list, unmapped_gene_list = self.fix_gene_names(gene_set)

        # Fetching genes
        if(gene_set): query_dictionary = {self.GeneField.FEATURE_TYPE:"gene", self.GeneField.GENE_ID:mapped_gene_list}
        else: query_dictionary = {self.GeneField.FEATURE_TYPE:"gene"}
        query_annset = self.get(query_dictionary)

        # Creating GenomicRegionSet
        result_grs = GenomicRegionSet("genes")
        for e in query_annset.gene_list:
            gr = e[self.GeneField.GENOMIC_REGION]
            gr.name = e[self.GeneField.GENE_ID]
            result_grs.add(gr)
        result_grs.merge()
        if(gene_set): return result_grs, unmapped_gene_list
        else: return result_grs

    def get_transcripts(self, gene_set = None):
        """Gets transcripts of genes. It returns a GenomicRegionSet with such transcripts. The id of each gene will be put in the NAME field of each GenomicRegion.

        *Keyword arguments:*

            - gene_set -- A set of genes to narrow the search.

        *Return:*

            - result_grs -- A GenomicRegionSet containing the exons.
            - unmapped_gene_list -- A list of genes that could not be mapped to an ENSEMBL ID.
        """

        # Fetching gene names
        mapped_gene_list = None
        unmapped_gene_list = None
        if gene_set: mapped_gene_list, unmapped_gene_list = self.fix_gene_names(gene_set)

        # Fetching exons
        if gene_set: query_dictionary = {self.GeneField.FEATURE_TYPE:"exon", self.GeneField.GENE_ID:mapped_gene_list}
        else: query_dictionary = {self.GeneField.FEATURE_TYPE:"exon"}
        query_annset = self.get(query_dictionary)

        # Creating GenomicRegionSet
        result_grs = GenomicRegionSet("exon")
        for e in query_annset.gene_list:
            gr = e[self.GeneField.GENOMIC_REGION]
            gr.name = e[self.GeneField.TRANSCRIPT_ID]
            result_grs.add(gr)

        if gene_set: return result_grs, unmapped_gene_list
        else: return result_grs

    def merge_rbs(self, rm_duplicate=False, asgene_organism=None, region_set=None, cutoff=0):
        """Merge the RNA binding regions which have overlap to each other and 
           combine their corresponding DNA binding regions.
        
        extend -> Define the extending length in basepair of each RNA binding regions
        perfect_match -> Merge only the exactly same RNA binding regions
        """
        # Merge RBS
        rna_merged = self.get_rbs()
        rna_merged.merge()
        # A dict: RBS as key, and GenomicRegionSet as its value
        new_dict = OrderedDict()

        for rbsm in rna_merged:
            regions = GenomicRegionSet(rbsm.toString())
            
            for rd in self:
                if rbsm.overlap(rd.rna):
                    regions.add(rd.dna)
            if rm_duplicate: 
                regions.remove_duplicates()
            if len(regions) > cutoff:
                new_dict[rbsm] = regions
                if asgene_organism:
                    try:
                        new_dict[rbsm] = new_dict[rbsm].gene_association(organism=asgene_organism)
                    except:
                        pass
                if region_set:
                    new_dict[rbsm].replace_region_name(regions=region_set)
            else: continue

        self.merged_dict = new_dict

 def region_sets(self,listA,listB):
     """ Setting two GenomicRegionSets as self.setA and self.setB for each case test. """
     self.setA = GenomicRegionSet('for Unit Test')
     for i in range(len(listA)):
         self.setA.add(GenomicRegion(chrom=listA[i][0], initial=listA[i][1], final=listA[i][2]))
     
     self.setB = GenomicRegionSet('for Unit Test')
     for i in range(len(listB)):
         self.setB.add(GenomicRegion(chrom=listB[i][0], initial=listB[i][1], final=listB[i][2]))

    def load_objects(self, is_bedgraph, verbose=False, test=False):
        """Load files and initialize object.

        *Keyword arguments:*

            - is_bedgraph -- Whether regions are in bedgraph format (default = False).
            - verbose -- Verbose output (default = False).
            - test -- Fetch only 10 regions form each BED files for test.
        """
        for i, t in enumerate(self.types):
            if verbose: print("Loading file ", self.files[self.names[i]], file = sys.stderr)
            
            if t not in ["regions", "genes"] and verbose:
                print("Cannot load objects", file=sys.stderr)
            
            if t == "regions":
                regions = GenomicRegionSet(self.names[i])
                if is_bedgraph:
                    regions.read_bedgraph(os.path.abspath(self.files[self.names[i]]))
                    
                else:
                    if test:
                        g = GenomicRegionSet(self.names[i])
                        g.read_bed(os.path.abspath(self.files[self.names[i]]))
                        regions.sequences = g.sequences[0:11]
                    else:
                        regions.read_bed(os.path.abspath(self.files[self.names[i]]))  # Here change the relative path into absolute path
                self.objectsDict[self.names[i]] = regions
            
            elif t == "genes":
                genes = GeneSet(self.names[i])
                genes.read(os.path.abspath(self.files[self.names[i]]))  # Here change the relative path into absolute path
                self.objectsDict[self.names[i]] = genes

    def get_promoters(self, promoterLength=1000, gene_set=None, unmaplist=False, variants=False):
        """
        Gets promoters of genes given a specific promoter length. It returns a GenomicRegionSet with such promoters.
        The ID of each gene will be put in the NAME field of each GenomicRegion.
        Each promoter includes also the coordinate of the 5' base pair, therefore each promoter actual
        length is promoterLength+1.

        *Keyword arguments:*

            - promoterLength -- The length of the promoter region.
            - gene_set -- A set of genes to narrow the search.
            - unmaplist -- If True than also return the unmappable genes list (default = False).

        *Return:*

            - result_grs -- A GenomicRegionSet containing the promoters.
            - unmapped_gene_list -- A list of genes that could not be mapped to an ENSEMBL ID.
        """

        # Fetching gene names
        mapped_gene_list = None
        unmapped_gene_list = None


        if gene_set: mapped_gene_list, unmapped_gene_list = self.fix_gene_names(gene_set)

        # Fetching genes

        if not variants: target = "gene"
        else: target = "transcript"
        if(gene_set): query_dictionary = {self.GeneField.FEATURE_TYPE:target, self.GeneField.GENE_ID:mapped_gene_list}
        else: query_dictionary = {self.GeneField.FEATURE_TYPE:target}
        
        query_annset = self.get(query_dictionary)

        # Creating GenomicRegionSet
        result_grs = GenomicRegionSet("promoters")
        for e in query_annset.gene_list:
            gr = e[self.GeneField.GENOMIC_REGION]
            if gr.orientation == "+":
                gr.final = gr.initial + 1
                gr.initial = gr.initial - promoterLength
            else:
                gr.initial = gr.final - 1
                gr.final = gr.initial + promoterLength + 1

            gr.name = e[self.GeneField.GENE_ID]
            result_grs.add(gr)
        if unmaplist: return result_grs, unmapped_gene_list
        else: return result_grs

 def get_dbs(self, sort=False, orientation=None, rm_duplicate=False):
     """Return GenomicRegionSet which contains all DNA binding sites"""
     dna_set = GenomicRegionSet(name="DNA_binding_sites")
     for rd in self.sequences:
         if not orientation:
             dna_set.add(rd.dna)
         else:
             if orientation == rd.orient:
                 dna_set.add(rd.dna)
             else: pass
     if sort: dna_set.sort()
     if rm_duplicate: dna_set.remove_duplicates()
     return dna_set

    def match_ms_tags(self,field):
        """Add more entries to match the missing tags of the given field. For example, there are tags for cell like 'cell_A' and 'cell_B' for reads, but no these tag for regions. Then the regions are repeated for each tags from reads to match all reads.

        *Keyword arguments:*

            - field -- Field to add extra entries.
        """
        
        # print(field)
        # print(self.fieldsDict)
        # check regions or reads have empty tag
        altypes = self.fieldsDict[field].keys()
        if "ALL" in altypes:
            altypes.remove("ALL")
            for name in self.fieldsDict[field]["ALL"]:
                # print(name)
                i = self.names.index(name)
                for t in altypes:
                    # print("\t"+t)
                    n = name+"_"+t
                    # print("\t\t"+n)
                    self.names.append(n)
                    self.types.append(self.types[i])
                    self.files[n] = self.files[name]
                    # types = self.get_types(name,skip_all=True)
                    # print("************")
                    # print(types)

                    for f in self.fields[3:]:
                        if f == field: 
                            try: self.fieldsDict[f][t].append(n)
                            except: self.fieldsDict[f][t] = [n]
                        else:
                            try: self.fieldsDict[f][self.get_type(name=name,field=f)].append(n)
                            except: self.fieldsDict[f][self.get_type(name=name,field=f)] = [n]
                    # for f in self.fieldsDict.keys():
                    #     for ty in types:
                    #         try: self.fieldsDict[f][ty].append(n)
                    #         except: pass
                    if self.types[i] == "regions":
                        g = GenomicRegionSet(n)
                        g.read_bed(self.files[name])
                        self.objectsDict[n] = g
                    self.trash.append(name)

def mode_3(exp_matrix, thresh, type_file):
    #remember value of bedgraph, ugly way
    score = {}
    for regions in exp_matrix.get_regionsets():
        for region in regions:
            if type_file=="ODIN":
              aux=(region.data).split("\t")
              aux=aux[-1].split(";")
              score[(region.chrom + ':' + str(region.initial) + '-' + str(region.final))] = float(region.data[-1])
            if type_file=="THOR":
              aux=(region.data).split(";")
              score[(region.chrom + ':' + str(region.initial) + '-' + str(region.final))] = float(aux[-1])
            else:
               score[(region.chrom + ':' + str(region.initial) + '-' + str(region.final))] = region.data
    
    
    for i, region in enumerate(exp_matrix.get_regionsets()):
        f = open("region_" + str(region.name) + ".data", 'w')
        region_set = GenomicRegionSet("")
        _, _, mappedGenes, _, gene_peaks_mapping = region_set.filter_by_gene_association_old(region.fileName, None, gene_file, genome_file, threshDist=thresh)

        avg_score = {} #score per peak
        genes = {}
        
        print('Consider row %s of exp. matrix, number of mapped genes is %s' %(i, mappedGenes), file=sys.stderr)
        for peak, gene_list in gene_peaks_mapping.items():            
            for gen in gene_list: #reverse mapping peak -> gene to gene -> peak
                if not gen:
                    continue
                genes[gen] = genes.get(gen, set())
                genes[gen].add(peak)
                
                avg_score[gen] = avg_score.get(gen, [])
                avg_score[gen].append(score[peak]) #join all scores of peaks assigned to a gen
        
        for gen in genes.keys():
            if options.metric == 'mean':
                avg = np.mean(avg_score[gen])
            elif options.metric == 'max':
                avg = np.max(avg_score[gen])
            print(gen, avg, ", ".join(str(t) for t in genes[gen]), sep='\t', file = f)
        
        f.close()

    def get_tts(self, gene_set=None):
        """Gets TTS(Transcription termination site) of genes. It returns a GenomicRegionSet with such TTS. The ID of each gene will be put in the NAME field of each GenomicRegion.

        *Keyword arguments:*

            - gene_set -- A set of genes to narrow the search.
        
        *Return:*

            - result_grs -- A GenomicRegionSet containing TTS.
            - unmapped_gene_list -- A list of genes that could not be mapped to an ENSEMBL ID.
        """

        # Fetching gene names
        mapped_gene_list = None
        unmapped_gene_list = None
        if gene_set: mapped_gene_list, unmapped_gene_list = self.fix_gene_names(gene_set)

        # Fetching genes
        if gene_set:
            query_dictionary = {self.GeneField.FEATURE_TYPE: "gene", self.GeneField.GENE_ID: mapped_gene_list}
        else:
            query_dictionary = {self.GeneField.FEATURE_TYPE: "gene"}
        query_annset = self.get(query_dictionary)

        # Creating GenomicRegionSet
        result_grs = GenomicRegionSet("TTS")
        for e in query_annset.gene_list:
            gr = e[self.GeneField.GENOMIC_REGION]
            if gr.orientation == "+":
                gr.initial = gr.initial
                gr.final = gr.initial + 1
            else:
                gr.initial = gr.final - 1
                gr.final = gr.final
            gr.name = e[self.GeneField.GENE_ID]
            result_grs.add(gr)
        result_grs.merge()
        if gene_set:
            return result_grs, unmapped_gene_list
        else:
            return result_grs

    def get_exons(self, start_site=False, end_site=False, gene_set=None, merge=True):
        """Gets exons of genes. It returns a GenomicRegionSet with such exons. The id of each gene will be put in the NAME field of each GenomicRegion.

        *Keyword arguments:*

            - start_site -- Whether to relocate the start sites.
            - end_site -- Whether to relocate the end sites.
            - gene_set -- A set of genes to narrow the search.

        *Return:*

            - result_grs -- A GenomicRegionSet containing the exons.
            - unmapped_gene_list -- A list of genes that could not be mapped to an ENSEMBL ID.
        """

        # Fetching gene names
        mapped_gene_list = None
        unmapped_gene_list = None
        if gene_set: mapped_gene_list, unmapped_gene_list = self.fix_gene_names(gene_set)

        # Fetching exons
        if gene_set:
            query_dictionary = {self.GeneField.FEATURE_TYPE: "exon", self.GeneField.GENE_ID: mapped_gene_list}
        else:
            query_dictionary = {self.GeneField.FEATURE_TYPE: "exon"}
        query_annset = self.get(query_dictionary)

        # Creating GenomicRegionSet
        result_grs = GenomicRegionSet("exon")
        for e in query_annset.gene_list:
            gr = e[self.GeneField.GENOMIC_REGION]
            # gr.name = e[self.GeneField.GENE_ID]
            gr.name = e[self.GeneField.TRANSCRIPT_ID]
            result_grs.add(gr)
        if start_site:
            result_grs.relocate_regions("leftend", left_length=1, right_length=1)
        elif end_site:
            result_grs.relocate_regions("rightend", left_length=1, right_length=1)
        if merge: result_grs.merge()
        if gene_set:
            return result_grs, unmapped_gene_list
        else:
            return result_grs

 def test_get_genome_data(self):
     """hg19"""
     result = GenomicRegionSet("hg19")
     result.get_genome_data(organism="hg19")
     self.assertEqual(len(result), 23)
     """hg19, with Mitochondria chromosome"""
     result = GenomicRegionSet("hg19")
     result.get_genome_data(organism="hg19",chrom_M=True)
     self.assertEqual(len(result), 24)

    def get_dbs(self, sort=False, orientation=None, rm_duplicate=False, dbd_tag=False):
        """Return GenomicRegionSet which contains all DNA binding sites"""
        dna_set = GenomicRegionSet(name="DNA_binding_sites")
        if len(self) == 0: return dna_set
        for rd in self.sequences:
            if dbd_tag:
                dbs = GenomicRegion(chrom=rd.dna.chrom, initial=rd.dna.initial, final=rd.dna.final,
                                    name=rd.rna.str_rna(), orientation=rd.dna.orientation, 
                                    data=rd.score)
            else:
                dbs = GenomicRegion(chrom=rd.dna.chrom, initial=rd.dna.initial, final=rd.dna.final,
                                    name=rd.dna.name, orientation=rd.dna.orientation, 
                                    data=rd.score)

            if not orientation:
                dna_set.add(dbs)
            else:
                if orientation == rd.orient:
                    dna_set.add(dbs)
                else: pass
        if sort: dna_set.sort()
        if rm_duplicate: dna_set.remove_duplicates()
        return dna_set

 def test_filter_tts(self):
 	txp = RNADNAInteractionSet(organism="hg19", filename=sample_txp)
 	g = GenomicRegionSet("g")
 	s = GenomicRegion(chrom="chr2", initial=74000000, final=75000000)
 	g.add(s)
 	result = txp.count_tts(g)

class TestGenomicRegionSet(unittest.TestCase):
    
    def region_sets(self,listA,listB):
        """ Setting two GenomicRegionSets as self.setA and self.setB for each case test. """
        self.setA = GenomicRegionSet('for Unit Test')
        for i in range(len(listA)):
            self.setA.add(GenomicRegion(chrom=listA[i][0], initial=listA[i][1], final=listA[i][2]))
        
        self.setB = GenomicRegionSet('for Unit Test')
        for i in range(len(listB)):
            self.setB.add(GenomicRegion(chrom=listB[i][0], initial=listB[i][1], final=listB[i][2]))
    
    def test_extend(self):
        """
        Two empty sets
        A : none 
        R : none
        """
        self.region_sets([],
                         [])
        self.setA.extend(100,100)
        self.assertEqual(len(self.setA.sequences), 0)
        """
        One region
        A :   -----
        R : ---------
        """
        self.region_sets([['chr1',5,10]],
                         [])
        result = self.setA
        result.extend(4,4)
        self.assertEqual(len(result), 1)
        self.assertEqual(result[0].initial, 1)
        self.assertEqual(result[0].final, 14)
        """
        Many region
        A :   -----   ------         -----    -----
        R : --------=---------     ------------------
        """
        self.region_sets([['chr1',5,10],['chr1',15,20],['chr1',40,50],['chr1',65,75]],
                         [])
        result = self.setA
        result.extend(5,5)
        self.assertEqual(len(result), 4)
        self.assertEqual(result[0].initial, 0)
        self.assertEqual(result[0].final, 15)
        self.assertEqual(result[1].initial, 10)
        self.assertEqual(result[1].final, 25)
        self.assertEqual(result[2].initial, 35)
        self.assertEqual(result[2].final, 55)
        self.assertEqual(result[3].initial, 60)
        self.assertEqual(result[3].final, 80)
        """
        Many region in different chromosome
        A :   -----   ------         -----    -----
        R : none
        """
        self.region_sets([['chr1',5,10],['chr2',15,20],['chr3',40,50],['chr4',65,75]],
                         [])
        result = self.setA
        result.extend(5,5)
        self.assertEqual(len(result), 4)
        self.assertEqual(result[0].initial, 0)
        self.assertEqual(result[0].final, 15)
        self.assertEqual(result[0].chrom, 'chr1')
        self.assertEqual(result[1].initial, 10)
        self.assertEqual(result[1].final, 25)
        self.assertEqual(result[1].chrom, 'chr2')
        self.assertEqual(result[2].initial, 35)
        self.assertEqual(result[2].final, 55)
        self.assertEqual(result[2].chrom, 'chr3')
        self.assertEqual(result[3].initial, 60)
        self.assertEqual(result[3].final, 80)
        self.assertEqual(result[3].chrom, 'chr4')
        """
        One region
        A :   -----
        R : ---------
        """
        self.region_sets([['chr1',100,200]],
                         [])
        result = self.setA
        result.extend(10,10,percentage=True)
        self.assertEqual(len(result), 1)
        self.assertEqual(result[0].initial, 90)
        self.assertEqual(result[0].final, 210)
        
    def test_sort(self):
        self.region_sets([['chr1',15,20],['chr1',40,50],['chr1',65,75],['chr1',5,10]],
                         [])
        self.setA.sort()
    
    def test_intersect(self):
        """
        Two empty sets
        A : none 
        B : none
        R : none
        """
        self.region_sets([],
#.........这里部分代码省略.........

from rgt.GenomicRegionSet import *
from rgt.ExperimentalMatrix import *
#from fisher import pvalue
import scipy.stats


outdir=""

back=False
designFile = sys.argv[1]
anotationPath = sys.argv[2]
randomize = int(sys.argv[3])
backGroundPeaks=False
if len(sys.argv) > 4:
  backGroundPeaksName = sys.argv[4]
  backBed=GenomicRegionSet("BACK")
  backBed.read_bed(backGroundPeaksName)  
  backGroundPeaks=True
   
distance=50000
if len(sys.argv) > 5:
  distance=len(sys.argv[5])

if len(sys.argv) > 6:
  outdir=sys.argv[6]


genomeFile=anotationPath+"chrom.sizes"
geneFile=anotationPath+"association_file.bed"

exps=ExperimentalMatrix()

def load_exon_sequence(bed, directory, genome_path):
    """Load the exon sequence from the the transcripts. 
    Input BED format should contain:
        blockCount - The number of blocks (exons) in the BED line.
        blockSizes - A comma-separated list of the block sizes.
        blockStarts - A comma-separated list of block starts. 
        see details: http://genome.ucsc.edu/FAQ/FAQformat#format1

    Output:
        Each FASTA file represants a transcript and contains all the exons within the file.

    """
    regionset = GenomicRegionSet("bed")
    regionset.read_bed(bed)
    regionset.sort()

    
    genome = pysam.Fastafile(genome_path)
    
    try:
        if len(regionset.sequences[0].data.split("\t")) == 7: 
            blockinfor = True
            no_exon = False
    except:
        blockinfor = False
        regionset.sequences.sort(key=lambda g: g.name)
        no_exon = True

    if blockinfor:
        
        for gr in regionset:
            if not gr.name:
                print("Error: For fetching exon sequences, please define the transcript name.")
                sys.exit()
            else:
                if not os.path.exists(directory):
                    os.makedirs(directory)
                f = open(os.path.join(directory, gr.name+".fa"), "w")
                data = gr.data.split("\t")
                #print(len(data))
                if len(data) == 7:
                    #print(data)
                    n = int(data[4])
                    
                    blocks = [ int(b) for b in filter(None, data[5].split(",")) ]
                    starts = [ int(s) for s in filter(None, data[6].split(",")) ]
                    printstr = []

                    for i in range(n):
                        start = gr.initial + starts[i]
                        end = start + blocks[i]
                        if no_exon and i == 0:
                            ex = ""
                        elif gr.orientation == "-":
                            ex = "exon:"+str(n-i)
                        else:
                            ex = "exon:"+str(i+1)

                        if gr.orientation == "-":
                            seq = Seq(genome.fetch(gr.chrom, start-1, end-1), IUPAC.unambiguous_dna)
                            seq = seq.reverse_complement()
                            p = [ ">"+ " ".join([ gr.name, 
                                                  ex, 
                                                  "_".join(["REGION",gr.chrom,
                                                            str(start),str(end), 
                                                            gr.orientation]) ]),
                                  seq ]
                            
                            printstr.append(p)
                            

                        else:
                            p = [ ">"+ " ".join([gr.name, ex, 
                                  "_".join(["REGION",gr.chrom,str(start),str(end), gr.orientation]) ]),
                                  genome.fetch(gr.chrom, start-1, end-1)
                                ]
                            printstr.append(p)
                            

                    if gr.orientation == "-": printstr = printstr[::-1]
                    for i in range(n):
                        print(printstr[i][0], file=f)
                        print(printstr[i][1], file=f)
                        

                else:
                    print("Warning: The given regions have no block information, please try write_bed_blocks")
                f.close()
    else:
        pre_id = ""
        for gr in regionset:
            if not gr.name: 
                gr.name = gr.toString()

            if pre_id == "": 
                pre_id = gr.name
                z = GenomicRegionSet(gr.name)
                z.add(gr)
            elif gr.name == pre_id:
                z.add(gr)
#.........这里部分代码省略.........