How to use iEnd method in fMBT
Best Python code snippet using fMBT_python
alternative_splicing.py
Source:alternative_splicing.py 
1import sys2import csv3from collections import defaultdict4from decimal import *5#Esta version tiene una modificacion para que la variantes de splicing que ocurran en el ultiumo exon sean bien clasificadas.6def main(final_table, genecode_bed):7 8 getcontext().prec = 209 10 istart_coverage = defaultdict(int)11 iend_coverage = defaultdict(int) 12 alt_start_end = defaultdict(list)13 14 intron_skipped_exons = {}15 intron_retention_exons = {}16 non_canonical = []17 18 TOTAL_intron = defaultdict(list)19 TOTAL_intron_set = set([])20 21 exons = set([])22 exons_coverage = defaultdict(list)23 24 introns_genecode = [] 25 26 27 for row in csv.reader(open(genecode_bed), delimiter = '\t'):28 qstarts = map (int, row[11].strip(",").split(",")) 29 blocksizes = map(int, row[10].strip(",").split(","))30 start = int(row[1])31 strand = row[5]32 bn = int(row[9])33 chr = row[0]34 35 blocks_number = len(blocksizes)36 blockcount = 037 for q, b in zip(qstarts, blocksizes):38 exon_pos = "middle"39 blockcount += 140 if blockcount == 1:41 exon_pos = "first"42 if blockcount == blocks_number:43 exon_pos = "last"44 45 estart = start+q46 eend = start+q+b47 exon = chr + ":" + str(estart) + strand + str(eend)48 exons.add((exon, chr, strand, estart, eend, exon_pos))49 50 for row in csv.reader(open(final_table), delimiter = ' '):51 intron = row[0]52 chr = row[1]53 strand = row[2]54 istart = int(row[3])55 iend = int(row[4])56 ilength = row[5]57 dn = row[6]58 dn_type = row[7]59 dn_type_score = row[8]60 bodymap_coverage = int(row[9])61 gm12878_coverage = int(row[10])62 hg19_cDNA_coverage = int(row[11])63 hg19_EST_coverage = int(row[12])64 mm9_cDNA_coverage = int(row[13])65 mm9_EST_coverage = int(row[14])66 genecode_coverage = int(row[15])67 tissues_coverage = int(row[16])68 n_tissues = row[17]69 tissues = row[18]70 71 istart_coverage[chr+strand+str(istart)] += bodymap_coverage + gm12878_coverage + tissues_coverage72 iend_coverage[chr+strand+str(iend)] += bodymap_coverage + gm12878_coverage + tissues_coverage73 74 alt_start_end[chr+strand+str(istart)+"s"].append(row)75 alt_start_end[chr+strand+str(iend)+"e"].append(row)76 77 78# if dn !="GTAG" and dn != "GCAG" and dn != "ATAC":79 if 0==0:80 81 non_canonical.append(row)82# else:83 TOTAL_intron[chr+strand].append(row)84 TOTAL_intron_set.add(intron)85 86 87 for e in exons:88 exon = e[0]89 chr = e[1]90 strand = e[2]91 estart = e[3]92 eend = e[4]93 exon_pos = e[5] 94 coverages = []95 coverages.append(iend_coverage[chr+strand+str(estart)])96 coverages.append(istart_coverage[chr+strand+str(eend)])97 98 if 0 in coverages:99 if exon_pos=="first" or exon_pos=="last": 100 coverages.remove(0)101 else:102 coverages = [0]103 exon_coverage = float(sum(coverages))/float(len(coverages))104 exons_coverage[chr+strand].append((exon, chr, strand, estart, eend, exon_pos, exon_coverage))105 106 107 for row in non_canonical:108 109 intron = row[0]110 chr = row[1]111 strand = row[2]112 istart = int(row[3])113 iend = int(row[4])114 ilength = row[5]115 dn = row[6]116 dn_type = row[7]117 dn_type_score = row[8]118 bodymap_coverage = int(row[9])119 gm12878_coverage = int(row[10])120 hg19_cDNA_coverage = int(row[11])121 hg19_EST_coverage = int(row[12])122 mm9_cDNA_coverage = int(row[13])123 mm9_EST_coverage = int(row[14])124 genecode_coverage = int(row[15])125 tissues_coverage = int(row[16])126 n_tissues = row[17]127 tissues = row[18]128 129 skipped_exons = []130 retention_exons = []131 for e in exons_coverage[chr+strand]:132 exon = e[0]133 chr = e[1]134 strand = e[2]135 estart = e[3]136 eend = e[4]137 exon_pos = e[5] 138 exon_coverage = e[6]139# if exon_pos == "middle": #Aqui me equiboque!!, los last tambien es mejor considerarlos como exon skiped!140 if istart<=estart and eend<=iend:141 skipped_exons.append((chr, strand, estart, eend, str(round(Decimal(exon_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) ))142 if estart<=istart and iend<=eend:143 retention_exons.append((chr, strand, estart, eend))144 145 intron_retention_exons[intron] = retention_exons 146 intron_skipped_exons[intron] = skipped_exons147 148 149 for row in non_canonical:150 151 152 intron = row[0]153 chr = row[1]154 strand = row[2]155 istart = int(row[3])156 iend = int(row[4])157 ilength = row[5]158 dn = row[6]159 dn_type = row[7]160 dn_type_score = row[8]161 bodymap_coverage = int(row[9])162 gm12878_coverage = int(row[10])163 hg19_cDNA_coverage = int(row[11])164 hg19_EST_coverage = int(row[12])165 mm9_cDNA_coverage = int(row[13])166 mm9_EST_coverage = int(row[14])167 genecode_coverage = int(row[15])168 tissues_coverage = int(row[16])169 n_tissues = row[17]170 tissues = row[18]171 172 info_start = []173 174 if alt_start_end.has_key(chr+strand+str(istart)+"s"):175 info_start = alt_start_end[chr+strand+str(istart)+"s"]176 177 tuple_info_start = map(lambda x: tuple(x), info_start)178 179 set_info_start = set(tuple_info_start)180 181 info_end = []182 183 184 if alt_start_end.has_key(chr+strand+str(iend)+"e"):185 info_end = alt_start_end[chr+strand+str(iend)+"e"]186 187 tuple_info_end = map(lambda x: tuple(x), info_end)188 189 set_info_end = set(tuple_info_end)190 191 alt_info = (set_info_start | set_info_end) - set([tuple(row)])192 193 alt_introns = []194 alt_skipper_introns = []195 alt_no_skipper_introns = []196 alt_exon_variant_introns = []197 198 199 alt_canonical = False #indicador para saber si hay un intron canonico con el donor o aceptor del no-canonico200 201 for alt in alt_info:202 alt_intron = alt[0]203 alt_chr = alt[1]204 alt_strand = alt[2]205 alt_istart = int(alt[3])206 alt_iend = int(alt[4])207 alt_dn = alt[6]208 alt_bodymap_coverage = int(alt[9])209 alt_gm12878_coverage = int(alt[10])210 alt_tissues_coverage = int(alt[16]) 211 skipped_exons = intron_skipped_exons[intron]212 213 if alt_dn == "GTAG" or alt_dn == "GCAG" or alt_dn == "ATAC" :214 alt_canonical = True215 216 alt_skipped_exons = []217 218 for e in exons_coverage[chr+strand]:219 exon = e[0]220 chr = e[1]221 strand = e[2]222 estart = e[3]223 eend = e[4]224 exon_pos = e[5] 225 exon_coverage = e[6]226# if exon_pos == "middle":227 if alt_istart<=estart and alt_iend>=eend:228 alt_skipped_exons.append((chr, strand, estart, eend, str(round(Decimal(exon_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) ))229 is_alt_no_skipper_intron = False 230 for e in skipped_exons: #Modulo para encontrar las variantes de intrones que son resultado no saltarse un intron231 if (alt_istart in e) or (alt_iend in e):232 is_alt_no_skipper_intron = True233 if is_alt_no_skipper_intron == True and (alt_istart != istart or alt_iend != iend) :234 if strand == "+":235 alt_no_skipper_introns.append(alt_intron + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(alt_istart - istart) + "|" + str(alt_iend - iend))236 if strand == "-":237 alt_no_skipper_introns.append(str(alt_intron) + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(iend - alt_iend) + "|" + str(istart - alt_istart)) 238 if set(alt_skipped_exons) == set(skipped_exons) and (alt_istart != istart or alt_iend != iend): #Modulo para en contrar variantes de intrones que son producidas por saltarse mas exones que el no canonico239 if strand == "+":240 alt_exon_variant_introns.append(alt_intron + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(alt_istart - istart) + "|" + str(alt_iend - iend))241 if strand == "-":242 alt_exon_variant_introns.append(str(alt_intron) + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(iend - alt_iend) + "|" + str(istart - alt_istart)) 243 244 if len(set(alt_skipped_exons) -set(skipped_exons))!=0 and (alt_istart != istart or alt_iend != iend): #Modulo para en contrar variantes de intrones que son producidas por saltarse mas exones que el no canonico245 if strand == "+":246 alt_skipper_introns.append(alt_intron + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(alt_istart - istart) + "|" + str(alt_iend - iend))247 if strand == "-":248 alt_skipper_introns.append(str(alt_intron) + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(iend - alt_iend) + "|" + str(istart - alt_istart)) 249 250 if strand == "+" and (alt_istart != istart or alt_iend != iend): #Modulo para encontrar el total de las variates 5' y 3'251 alt_introns.append(alt_intron + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(alt_istart - istart) + "|" + str(alt_iend - iend))252 if strand == "-" and (alt_istart != istart or alt_iend != iend):253 alt_introns.append(str(alt_intron) + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(iend - alt_iend) + "|" + str(istart - alt_istart))254 shift = []255 non_canonical_shift = []256 retention_exons = []257 retention_exons_names = []258 skipped_exons = []259 skipped_exons_names = []260 261 try:262 retention_exons = intron_retention_exons[intron]263 except KeyError:264 pass 265 266 for e in retention_exons:267 chr = e[0]268 strand = e[1]269 estart = e[2]270 eend = e[3]271 name = chr + ":" + str(estart) + strand + str(eend)272 retention_exons_names.append(name)273 274 try:275 skipped_exons = intron_skipped_exons[intron]276 except KeyError:277 pass 278 279 for e in skipped_exons:280 chr = e[0]281 strand = e[1]282 estart = e[2]283 eend = e[3]284 fold = e[4]285 name = chr + ":" + str(estart) + strand + str(eend)286 skipped_exons_names.append((name + "|" + fold)) 287 288 289 if alt_canonical == False: #solo busca si hay shifts cuando no se encontro un intron canonico alternativo290 TOTAL_intron_info = TOTAL_intron[chr+strand]291 for alt in TOTAL_intron_info:292 alt_intron = alt[0]293 alt_chr = alt[1]294 alt_strand = alt[2]295 alt_istart = int(alt[3])296 alt_iend = int(alt[4])297 alt_dn = alt[6]298 alt_bodymap_coverage = int(alt[9])299 alt_gm12878_coverage = int(alt[10])300 301 if (alt_istart-10) <= istart and istart <= (alt_istart+10) and (alt_iend-10) <= iend and iend <= (alt_iend + 10):302 if (alt_dn == "GTAG" or alt_dn == "GCAG" or alt_dn == "ATAC") and (alt_istart != istart or alt_iend != iend) : 303 if strand== "+":304 shift.append(alt_intron + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(alt_istart - istart) + "|" + str(alt_iend - iend))305 if strand== "-":306 shift.append(str(alt_intron) + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(iend - alt_iend) + "|" + str(istart - alt_istart))307 elif (alt_istart - istart) !=0 and (alt_iend - iend) != 0 and (alt_istart != istart or alt_iend != iend): #Para evitar que se reporte el intron no canonico como shift de si mismo, o que alt 3 o 5308 if strand== "+":309 non_canonical_shift.append(alt_intron + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(alt_istart - istart) + "|" + str(alt_iend - iend))310 if strand== "-":311 non_canonical_shift.append(str(alt_intron) + "|" + alt_dn + "|" + str(round(Decimal(alt_bodymap_coverage + alt_gm12878_coverage + alt_tissues_coverage)/Decimal(bodymap_coverage + gm12878_coverage + tissues_coverage), 3)) + "|" + str(iend - alt_iend) + "|" + str(istart - alt_istart))312 313 if bodymap_coverage + gm12878_coverage == 0:314 retention_exons_names = ""315 316 alt_table = [ ','.join(retention_exons_names), ','.join(skipped_exons_names), ','.join(alt_introns), ','.join(alt_no_skipper_introns), ','.join(alt_skipper_introns), ','.join(alt_exon_variant_introns), ','.join(shift), ','.join(non_canonical_shift)]317 318 print "\t".join(row) + "\t" + "\t".join(["NO" if x=="" else x for x in alt_table])319 320 321 322#python ~/my_src/Analisys/Alternative_splicing/alternative_splicing.py FINAL_TABLE ~/db/transcriptome/hg19/Gene_models/gencode/v11/gencode.v11.annotation.bed12323if __name__ == '__main__':...
BinWindUtils.py
Source:BinWindUtils.py
1import numpy as np2from astropy.io import ascii3from astropy.table import Table4#from Constants import Constants5import athena_read as ar6def get_slice_val(fn,level,slice_val=0.0):7 d = ar.athdf(fn,quantities=[],level=level,subsample=True)8 return d['x3v'][ np.argmin(np.abs(d['x3v'] - slice_val) ) ]9def read_data(fn,level,10 get_slice=True,slice_dir='z',slice_val=0.0,11 lim=48.):12 print("retrieving data with level =", level," within in limit=",lim)13 14 if(get_slice):15 sv = get_slice_val(fn,level,slice_val=slice_val)16 print("slicing at ",slice_dir,"=",sv)17 18 if(slice_dir=='z'):19 d = ar.athdf(fn,20 x1_min=-lim,x1_max=lim,21 x2_min=-lim,x2_max=lim,22 x3_min=sv,x3_max=sv,23 subsample=True,24 level=level)25 elif(slice_dir=='y'):26 d = ar.athdf(fn,27 x1_min=-lim,x1_max=lim,28 x2_min=sv,x2_max=sv,29 x3_min=-lim,x3_max=lim,30 subsample=True,31 level=level) 32 elif(slice_dir=='x'):33 d = ar.athdf(fn,34 x1_min=sv,x1_max=sv,35 x2_min=-lim,x2_max=lim,36 x3_min=-lim,x3_max=lim,37 subsample=True,38 level=level)39 40 else:41 d = ar.athdf(fn,42 x1_min=-lim,x1_max=lim,43 x2_min=-lim,x2_max=lim,44 x3_min=-lim,x3_max=lim,45 subsample=True,46 level=level)47 48 x1 = -0.549 x2 = 0.550 #dx = d['x1v'][1]-d['x1v'][0]51 #dA = dx**252 #dV = dx**353 54 d['x']=np.broadcast_to(d['x1v'],(len(d['x3v']),len(d['x2v']),len(d['x1v'])) )55 d['y']=np.swapaxes(np.broadcast_to(d['x2v'],(len(d['x3v']),len(d['x1v']),len(d['x2v'])) ),1,2)56 d['z']=np.swapaxes(np.broadcast_to(d['x3v'],(len(d['x1v']),len(d['x2v']),len(d['x3v'])) ) ,0,2 )57 58 d1 = np.sqrt((d['x']-x1)**2 +d['y']**2 + d['z']**2 )59 d2 = np.sqrt((d['x']-x2)**2 +d['y']**2 + d['z']**2 )60 d['PhiEff'] = -0.5/d1 -0.5/d2 - 0.5*(d['x']**2 + d['y']**2)61 62 # torque on binary due to wind63 d['torque_dens_z'] = x1*(0.5*d['rho']/d1**3 * d['y']) + x2*(0.5*d['rho']/d2**3 * d['y'])64 65 d['vx'] = d['vel1'] - 1.0*d['y'] # rot frame "+ Omega x r"66 d['vy'] = d['vel2'] + 1.0*d['x']67 d['vz'] = d['vel3']68 69 print("data has shape",d['rho'].shape)70 71 return d72# mdot = A * rho * v (cm^2 * g/cm^3 * cm/s)73def get_fluxes(d,offset=0):74 75 flux = {}76 77 istart = offset78 iend = len(d['x1v'])-1-offset79 80 print("box limits = ",d['x1v'][istart],d['x1v'][iend] )81 dx = d['x1v'][1]-d['x1v'][0]82 dA = dx**283 flux_mdot = np.zeros(6)84 flux_jdot_z = np.zeros(6)85 86 87 # -z88 flux_mdot[0] = np.sum( d['rho'][istart,istart:iend,istart:iend]*(-d['vz'][istart,istart:iend,istart:iend])*dA )89 flux_jdot_z[0] = np.sum( d['rho'][istart,istart:iend,istart:iend]*90 (d['x'][istart,istart:iend,istart:iend]*d['vy'][istart,istart:iend,istart:iend] 91 - d['y'][istart,istart:iend,istart:iend]*d['vx'][istart,istart:iend,istart:iend])*92 (-d['vz'][istart,istart:iend,istart:iend])*dA)93 # +z 94 flux_mdot[1] = np.sum( d['rho'][iend,istart:iend,istart:iend]*d['vz'][iend,istart:iend,istart:iend]*dA )95 flux_jdot_z[1] = np.sum( d['rho'][iend,istart:iend,istart:iend]*96 (d['x'][iend,istart:iend,istart:iend]*d['vy'][iend,istart:iend,istart:iend] 97 - d['y'][iend,istart:iend,istart:iend]*d['vx'][iend,istart:iend,istart:iend])*98 d['vz'][iend,istart:iend,istart:iend]*dA)99 100 # -y101 flux_mdot[2] = np.sum( d['rho'][istart:iend,istart,istart:iend]*(-d['vy'][istart:iend,istart,istart:iend])*dA )102 flux_jdot_z[2] = np.sum( d['rho'][istart:iend,istart,istart:iend]*103 (d['x'][istart:iend,istart,istart:iend]*d['vy'][istart:iend,istart,istart:iend] 104 - d['y'][istart:iend,istart,istart:iend]*d['vx'][istart:iend,istart,istart:iend])*105 (-d['vy'][istart:iend,istart,istart:iend])*dA)106 # +y107 flux_mdot[3] = np.sum( d['rho'][istart:iend,iend,istart:iend]*d['vy'][istart:iend,iend,istart:iend]*dA )108 flux_jdot_z[3] = np.sum( d['rho'][istart:iend,iend,istart:iend]*109 (d['x'][istart:iend,iend,istart:iend]*d['vy'][istart:iend,iend,istart:iend] 110 - d['y'][istart:iend,iend,istart:iend]*d['vx'][istart:iend,iend,istart:iend])*111 d['vy'][istart:iend,iend,istart:iend]*dA)112 113 # -x 114 flux_mdot[4] = np.sum( d['rho'][istart:iend,istart:iend,istart]*(-d['vx'][istart:iend,istart:iend,istart])*dA )115 flux_jdot_z[4] = np.sum( d['rho'][istart:iend,istart:iend,istart]*116 (d['x'][istart:iend,istart:iend,istart]*d['vy'][istart:iend,istart:iend,istart] 117 - d['y'][istart:iend,istart:iend,istart]*d['vx'][istart:iend,istart:iend,istart])118 *(-d['vx'][istart:iend,istart:iend,istart])*dA )119 # +x 120 flux_mdot[5] = np.sum( d['rho'][istart:iend,istart:iend,iend]*d['vx'][istart:iend,istart:iend,iend]*dA )121 flux_jdot_z[5] = np.sum( d['rho'][istart:iend,istart:iend,iend]*122 (d['x'][istart:iend,istart:iend,iend]*d['vy'][istart:iend,istart:iend,iend] 123 - d['y'][istart:iend,istart:iend,iend]*d['vx'][istart:iend,istart:iend,iend])124 *d['vx'][istart:iend,istart:iend,iend]*dA )125 flux['mdot'] = np.sum(flux_mdot)126 flux['mdot_fluxes'] = flux_mdot127 flux['Ldotz'] = np.sum(flux_jdot_z)128 flux['Ldotz_fluxes'] = flux_jdot_z129 130 return flux131def get_fluxes_sphere(d,radius):132 flux={}133 134 flux_mdot = np.zeros(6)135 flux_jdot_z = np.zeros(6)136 137 rin = np.where( np.sqrt(d['x']**2 + d['y']**2 + d['z']**2)<radius, True, False)138 cond_mz = ((rin[1:-1,1:-1,1:-1]==False) & (rin[2:,1:-1,1:-1]==True))139 cond_pz = ((rin[1:-1,1:-1,1:-1]==False) & (rin[0:-2,1:-1,1:-1]==True))140 cond_my = ((rin[1:-1,1:-1,1:-1]==False) & (rin[1:-1,2:,1:-1]==True))141 cond_py = ((rin[1:-1,1:-1,1:-1]==False) & (rin[1:-1,0:-2,1:-1]==True))142 cond_mx = ((rin[1:-1,1:-1,1:-1]==False) & (rin[1:-1,1:-1,2:]==True))143 cond_px = ((rin[1:-1,1:-1,1:-1]==False) & (rin[1:-1,1:-1,0:-2]==True))144 dx = d['x1v'][1]-d['x1v'][0]145 dA = dx**2146 d['l'] = d['x']*d['vy'] - d['y']*d['vx']147 flux_mdot[0] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_mz]*(-d['vz'][1:-1,1:-1,1:-1][cond_mz])*dA)148 flux_mdot[1] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_pz]*(d['vz'][1:-1,1:-1,1:-1][cond_pz])*dA)149 flux_mdot[2] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_my]*(-d['vy'][1:-1,1:-1,1:-1][cond_my])*dA)150 flux_mdot[3] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_py]*(d['vy'][1:-1,1:-1,1:-1][cond_py])*dA)151 flux_mdot[4] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_mx]*(-d['vx'][1:-1,1:-1,1:-1][cond_mx])*dA)152 flux_mdot[5] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_px]*(d['vx'][1:-1,1:-1,1:-1][cond_px])*dA)153 flux_jdot_z[0] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_mz]*(-d['vz'][1:-1,1:-1,1:-1][cond_mz])*dA154 *d['l'][1:-1,1:-1,1:-1][cond_mz])155 flux_jdot_z[1] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_pz]*(d['vz'][1:-1,1:-1,1:-1][cond_pz])*dA156 *d['l'][1:-1,1:-1,1:-1][cond_pz])157 flux_jdot_z[2] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_my]*(-d['vy'][1:-1,1:-1,1:-1][cond_my])*dA158 *d['l'][1:-1,1:-1,1:-1][cond_my])159 flux_jdot_z[3] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_py]*(d['vy'][1:-1,1:-1,1:-1][cond_py])*dA160 *d['l'][1:-1,1:-1,1:-1][cond_py])161 flux_jdot_z[4] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_mx]*(-d['vx'][1:-1,1:-1,1:-1][cond_mx])*dA162 *d['l'][1:-1,1:-1,1:-1][cond_mx])163 flux_jdot_z[5] = np.sum(d['rho'][1:-1,1:-1,1:-1][cond_px]*(d['vx'][1:-1,1:-1,1:-1][cond_px])*dA164 *d['l'][1:-1,1:-1,1:-1][cond_px])165 flux['mdot'] = np.sum(flux_mdot)166 flux['mdot_fluxes'] = flux_mdot167 flux['Ldotz'] = np.sum(flux_jdot_z)168 flux['Ldotz_fluxes'] = flux_jdot_z169 170 return flux171def sum_torque(d,fphi,offset=0):172 173 istart = offset174 iend = len(d['x1v'])-offset175 176 # sum torque = Ldot_grav177 dx = d['x1v'][1]-d['x1v'][0]178 dV = dx**3 179 cond = ((d['PhiEff']<-2*fphi) & (np.sqrt(d['x']**2 + d['y']**2 + d['z']**2)<1.2))[istart:iend,istart:iend,istart:iend]180 Ldot_grav = np.sum( np.ma.masked_where(cond,181 (d['torque_dens_z'][istart:iend,istart:iend,istart:iend])*dV) )182 return Ldot_grav183def sum_torque_sphere(d,fphi,radius):184 185 # sum torque = Ldot_grav186 dx = d['x1v'][1]-d['x1v'][0]187 dV = dx**3 188 cond = (((d['PhiEff']<-2*fphi) & (np.sqrt(d['x']**2 + d['y']**2 + d['z']**2)<1.2)) |189 (np.sqrt(d['x']**2 + d['y']**2 + d['z']**2)>radius))190 Ldot_grav = np.sum( np.ma.masked_where(cond,d['torque_dens_z']*dV) )191 return Ldot_grav192def get_mdot_gammas(d,fphi,offset=0):193 f = get_fluxes(d,offset)194 195 # l_loss 196 lbin = 0.25197 lloss = -f['Ldotz']/-f['mdot']198 gamma_loss = lloss/lbin 199 Ldot_grav = sum_torque(d,fphi,offset)200 201 # wind 202 # (Lgrav > 0 accelerates binary, decelerates wind)203 # (Lgrav < 0 decelarates binary, accelerates wind)204 lgrav = - Ldot_grav / f['mdot'] 205 gamma_grav = lgrav/lbin206 207 lwind = lloss - lgrav208 gamma_wind = lwind/lbin 209 210 return -f['mdot'], -f['Ldotz'], Ldot_grav, gamma_wind, gamma_grav, gamma_loss211def get_mdot_gammas_sphere(d,fphi,radius):212 f = get_fluxes_sphere(d,radius)213 214 # l_loss 215 lbin = 0.25216 lloss = -f['Ldotz']/-f['mdot']217 gamma_loss = lloss/lbin 218 Ldot_grav = sum_torque_sphere(d,fphi,radius)219 220 # wind 221 # (Lgrav > 0 accelerates binary, decelerates wind)222 # (Lgrav < 0 decelarates binary, accelerates wind)223 lgrav = - Ldot_grav / f['mdot'] 224 gamma_grav = lgrav/lbin225 226 lwind = lloss - lgrav227 gamma_wind = lwind/lbin 228 229 return -f['mdot'], -f['Ldotz'], Ldot_grav, gamma_wind, gamma_grav, gamma_loss230def get_vr10(d):231 dV=(d['x1v'][1]-d['x1v'][0])**3232 d['r'] = np.sqrt(d['x']**2 + d['y']**2 + d['z']**2)233 d['vr'] = d['vx']*d['x']/d['r'] + d['vy']*d['y']/d['r'] + + d['vz']*d['z']/d['r']234 sel10 = (d['r']>9.5) & (d['r']<10.5)235 vr10=np.sum(d['vr'][sel10]*d['rho'][sel10]*dV)/np.sum(d['rho'][sel10]*dV)...
ME_simulation.py
Source:ME_simulation.py
1import sys2import csv3from Bio import SeqIO4from Bio.Seq import Seq5from Bio.Alphabet import generic_dna6from random import randint, sample7from operator import itemgetter8from collections import defaultdict9Transcriptome = {}10Genome = {}11chr_sizes = {}12new_genome = {}13chrs = set([])14splice_sites = defaultdict(list)15def GenomeEditor(chr, micro_exon_istart, micro_exon_iend, strand, new_intron5, new_intron3, genome=Genome): #, splice_sites=splice_sites):16 if strand == "+":17 new_chr = Genome[chr][:micro_exon_iend-12] + new_intron3 + Genome[chr][micro_exon_iend:micro_exon_istart] + new_intron5 + Genome[chr][micro_exon_istart+6:]18 Genome[chr] = new_chr 19 if strand == "-":20 new_chr = Genome[chr][:micro_exon_istart-6] + new_intron5 + Genome[chr][micro_exon_istart:micro_exon_iend] + new_intron3 + Genome[chr][micro_exon_iend+12:]21 Genome[chr] = new_chr 22def Genomictabulator(fasta):23 24 print >> sys.stderr, "Cargando genoma en la memoria RAM ...", 25 f = open(fasta)26 for chrfa in SeqIO.parse(f, "fasta"):27 Genome[chrfa.id] = chrfa.seq28 29 print >> sys.stderr, "OK"30 f.close()31 32def Transcriptometabulator(genecode_fasta):33 34 print >> sys.stderr, "Cargando a fasta en la ram ...",35 36 for record in SeqIO.parse(genecode_fasta, "fasta"):37 id = str(record.id).split("|")[0]38 Transcriptome[id] = record.seq39 40 print >> sys.stderr, "OK"41 42 43 44def main(bed12):45 print >> sys.stderr, "Extrayendo intrones del bed12 ...",46 # row[0] chr47 # row[1] alignment start48 # row[2] alignment end49 # row[3] name 50 # row[4] 51 # row[5] strand52 # row[6] aligment start53 # row[7] aligment end54 # row[8] 55 # row[9] blocknum56 # row[10] blocksizes57 # row[11] qstarts 58 simulation_genome = open("simulation_genome.fa", 'w') 59 csv.field_size_limit(1000000000)60 n = 9961 min_intron_lenght = 80 62 introns = defaultdict(set)63 acepted_introns = set([])64 transcript_intron_info = {}65 for row in csv.reader(open(bed12), delimiter = '\t'): #Generando lista de intrones no solapantes66 67 68 qName = row[3]69 qstarts = map (int, row[11].strip(",").split(",")) 70 blocksizes = map(int, row[10].strip(",").split(","))71 start = int(row[1])72 strand = row[5]73 bn = int(row[9])74 chr = row[0]75 qstart = 076 for q1, q2, b in zip(qstarts, qstarts[1:], blocksizes):77 78 qstart = qstart + b79 tag_start = qstart - n80 tag_end = qstart + n81 istart = start + q1 + b82 iend = start + q283 ilen = iend - istart84 intron = row[0] + ":" + str(istart) + row[5] + str(iend) 85 intron = chr + ":" + str(istart) + strand + str(iend)86 ilength = iend - istart87 dn = Genome[chr][istart:(istart+2)] + Genome[chr][(iend-2):iend]88 intron5 = Genome[chr][istart:istart+6]89 intron3 = Genome[chr][iend-12:iend]90 91 if strand == "-":92 intron3 = Genome[chr][istart:istart+12]93 intron5 = Genome[chr][iend-6:iend]94 dn = dn.reverse_complement()95 dn = str(dn).upper()96 if ilen >= 250 and (dn=="GTAG" or dn=="GCAG" or dn=="ATAC"):97 introns[chr].add((intron, chr, istart, iend))98 splice_sites["intron5" + strand].append(intron5)99 splice_sites["intron3" + strand].append(intron3)100 for c in introns.items():101 chr, intron_list = c102 previous_end = 0103 for i in sorted(intron_list, key=itemgetter(2)):104 intron, chr, istart, iend = i105 if istart > previous_end:106 acepted_introns.add(intron)107 previous_end = iend108 for row in csv.reader(open(bed12), delimiter = '\t'):109 110 try:111 112 qName = row[3]113 seq = Transcriptome[qName]114 qstarts = map (int, row[11].strip(",").split(",")) 115 blocksizes = map(int, row[10].strip(",").split(","))116 start = int(row[1])117 strand = row[5]118 bn = int(row[9])119 chr = row[0]120 qstart = 0121 for q1, q2, b in zip(qstarts, qstarts[1:], blocksizes):122 123 qstart = qstart + b124 tag_start = qstart - n125 tag_end = qstart + n126 istart = start + q1 + b127 iend = start + q2128 ilen = iend - istart129 intron = row[0] + ":" + str(istart) + row[5] + str(iend) 130 intron = chr + ":" + str(istart) + strand + str(iend)131 ilength = iend - istart132 intron_seq = Genome[chr][istart:iend]133 134 if strand == '+' : #Para los que aliniean en la hebra +135 136 if tag_start<0: #Precausiones generar buenos tag del primer y ultimo tag137 tag_start = 0 138 if tag_end>len(seq):139 tag_end=len(seq)140 tag = seq[tag_start:tag_end]141 142 if strand == '-' :143 intron_seq = intron_seq.reverse_complement()144 145 if tag_end>len(seq): #Para los que alinian en la hebra - es todo al inverso146 tag_end=len(seq) 147 tag = seq[-tag_end:-tag_start]148 if tag_start<=0: 149 tag = seq[-tag_end:]150 151 if len(tag) == 2*n: #and (intron in acepted_introns):152 info = intron_seq, min_intron_lenght, tag, chr, istart, iend, strand153 transcript_intron_info[intron] = info154 except KeyError:155 pass156 for intron in acepted_introns:157 try:158 intron_seq, min_intron_lenght, tag, chr, istart, iend, strand = transcript_intron_info[intron]159 core_intron = intron_seq[min_intron_lenght:-min_intron_lenght]160 micro_exon_start = randint(0, len(core_intron))161 micro_exon_lenght = randint(1,25)162 micro_exon_tag = tag[:n] + core_intron[micro_exon_start:micro_exon_start + micro_exon_lenght] + tag[n:]163 micro_exon_iend = istart + min_intron_lenght + micro_exon_start164 micro_exon_istart = micro_exon_iend + micro_exon_lenght165 if strand == "-":166 micro_exon_istart = iend - (min_intron_lenght + micro_exon_start + micro_exon_lenght) #Los di vuelta167 micro_exon_iend = micro_exon_istart + micro_exon_lenght168 intron5 = Genome[chr][micro_exon_istart:micro_exon_istart+6]169 intron3 = Genome[chr][micro_exon_iend-12:micro_exon_iend]170 171 if strand == "-":172 intron3 = Genome[chr][micro_exon_iend:micro_exon_iend+12].reverse_complement() #cambiado173 intron5 = Genome[chr][micro_exon_istart-6:micro_exon_istart].reverse_complement()174 new_intron5 = sample(splice_sites["intron5" + strand],1)[0]175 new_intron3 = sample(splice_sites["intron3" + strand],1)[0]176 GenomeEditor(chr, micro_exon_istart, micro_exon_iend, strand, new_intron5, new_intron3)177 print intron, micro_exon_tag, core_intron[micro_exon_start:micro_exon_start + micro_exon_lenght], micro_exon_start, micro_exon_start + micro_exon_lenght, intron3, intron5, micro_exon_iend, micro_exon_istart, new_intron3, new_intron5178 except KeyError:179 pass180 for i in Genome.items():181 chr, seq = i182 seq = str(seq)183 simulation_genome.write(">" + chr + "\n")184 simulation_genome.write(seq + "\n") 185if __name__ == '__main__':186 Genomictabulator(sys.argv[1])187 Transcriptometabulator(sys.argv[2])188 main (sys.argv[3])189#python my_src/Micro-exons/micro_exon_simulation.py db/genomes/hg19.fa db/trasncriptomes/Gene_models/gencode/v19/gencode.v19.pc_transcripts.fa db/trasncriptomes/Gene_models/gencode/v19/gencode.v19.annotation.bed190# if strand == "+":191#192# for n, c in zip(new_intron5, range(micro_exon_istart, micro_exon_istart+6)):193# new_genome[chr + ":" + c] = n194#195# for n, c in zip(new_intron3, range(micro_exon_iend-12, micro_exon_iend)):196# new_genome[chr + ":" + c] = n197#198# if strand == "-":199#200# for n, c in zip(new_intron3, range(micro_exon_istart, micro_exon_istart+12)):201# new_genome[chr + ":" + c] = n202#203# for n, c in zip(new_intron5, range(micro_exon_iend-6, micro_exon_iend)):...
Automation Testing Tutorials
Learn to execute automation testing from scratch with LambdaTest Learning Hub. Right from setting up the prerequisites to run your first automation test, to following best practices and diving deeper into advanced test scenarios. LambdaTest Learning Hubs compile a list of step-by-step guides to help you be proficient with different test automation frameworks i.e. Selenium, Cypress, TestNG etc.
LambdaTest Learning Hubs:
- JUnit Tutorial
- TestNG Tutorial
- Webdriver Tutorial
- WebDriverIO Tutorial
- Protractor Tutorial
- Selenium 4 Tutorial
- Jenkins Tutorial
- NUnit Tutorial
- Jest Tutorial
- Playwright Tutorial
- Cypress Tutorial
- PyTest Tutorial
YouTube
You could also refer to video tutorials over LambdaTest YouTube channel to get step by step demonstration from industry experts.
Run fMBT automation tests on LambdaTest cloud grid
Perform automation testing on 3000+ real desktop and mobile devices online.
Try LambdaTest Now !!
Get 100 minutes of automation test minutes FREE!!
Was this article helpful?
