本文整理汇总了Python中rgt.GenomicRegionSet类的典型用法代码示例。如果您正苦于以下问题:Python GenomicRegionSet类的具体用法?Python GenomicRegionSet怎么用?Python GenomicRegionSet使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了GenomicRegionSet类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的Python代码示例。
示例1: mode_2
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()
开发者ID:Marvin84,项目名称:reg-gen,代码行数:27,代码来源:geneOntologyFromBed.py
示例2: get_biotypes
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
开发者ID:CostaLab,项目名称:reg-gen,代码行数:29,代码来源:AnnotationSet.py
示例3: mode_3
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()
开发者ID:Marvin84,项目名称:reg-gen,代码行数:31,代码来源:geneOntologyFromBed.py
示例4: get_genes
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
开发者ID:jovesus,项目名称:reg-gen,代码行数:33,代码来源:AnnotationSet.py
示例5: get_transcripts
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
开发者ID:eggduzao,项目名称:reg-gen,代码行数:32,代码来源:AnnotationSet.py
示例6: merge_rbs
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
开发者ID:jovesus,项目名称:reg-gen,代码行数:33,代码来源:RNADNABindingSet.py
示例7: region_sets
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]))
开发者ID:eggduzao,项目名称:reg-gen,代码行数:9,代码来源:test_GenomicRegionSet.py
示例8: load_objects
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
开发者ID:Marvin84,项目名称:reg-gen,代码行数:33,代码来源:ExperimentalMatrix.py
示例9: get_promoters
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
开发者ID:eggduzao,项目名称:reg-gen,代码行数:50,代码来源:AnnotationSet.py
示例10: get_dbs
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
开发者ID:jovesus,项目名称:reg-gen,代码行数:13,代码来源:RNADNABindingSet.py
示例11: match_ms_tags
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)
开发者ID:Marvin84,项目名称:reg-gen,代码行数:44,代码来源:ExperimentalMatrix.py
示例12: mode_3
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()
开发者ID:jovesus,项目名称:reg-gen,代码行数:43,代码来源:genesFromBed.py
示例13: get_tts
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
开发者ID:CostaLab,项目名称:reg-gen,代码行数:42,代码来源:AnnotationSet.py
示例14: get_exons
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
开发者ID:CostaLab,项目名称:reg-gen,代码行数:43,代码来源:AnnotationSet.py
示例15: test_get_genome_data
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)
开发者ID:eggduzao,项目名称:reg-gen,代码行数:9,代码来源:test_GenomicRegionSet.py
示例16: get_dbs
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
开发者ID:eggduzao,项目名称:reg-gen,代码行数:23,代码来源:RNADNABindingSet.py
示例17: test_filter_tts
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)
开发者ID:Marvin84,项目名称:reg-gen,代码行数:6,代码来源:test_RNADNAInteractionSet.py
示例18: TestGenomicRegionSet
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([],
#.........这里部分代码省略.........
开发者ID:eggduzao,项目名称:reg-gen,代码行数:101,代码来源:test_GenomicRegionSet.py
示例19: int
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()
开发者ID:Marvin84,项目名称:reg-gen,代码行数:31,代码来源:geneAssociationZscore.py
示例20: load_exon_sequence
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)
#.........这里部分代码省略.........
开发者ID:Marvin84,项目名称:reg-gen,代码行数:101,代码来源:bed2fasta.py
注:本文中的rgt.GenomicRegionSet类示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论