Add v0.1 script

Emulatior_v0.1.py

@@ -0,0 +1,129 @@
+
+# coding: utf-8
+
+# In[16]:
+
+import numpy as np
+import random
+import copy
+
+#N_observ=10
+#N_subclones=2
+N_observ=int(raw_input('Size of emulated data?'))
+N_subclones=int(raw_input('Number of emulated subclones?'))
+norm_data=np.zeros((N_observ,4))
+
+# DP Emulate
+for i in norm_data:
+    i[0]=1
+    i[3]=1
+    candidat=random.randint(1, N_observ*100)
+    while candidat in norm_data[:,1]:
+        candidat=random.randint(1, N_observ*100)
+    i[1]=candidat
+    i[2]=int(random.normalvariate(100, 20))
+norm_data=norm_data[np.argsort(norm_data[:,1])]
+
+# Subclones CNA Emulate
+SC=[]
+BSC=[]
+VSC=[]
+snv_dict={}
+change_CN=[0, 0.5, 1.5,2]
+for i in range(N_subclones):
+    SC.append(copy.deepcopy(norm_data))
+    BSC.append(copy.deepcopy(norm_data))
+    for j in BSC[i]:
+        j[3]=j[2]
+        j[2]=0
+        a=random.random()
+        if a<=0.5:
+            j[2]=int(random.normalvariate(j[3]/2, 2))
+    coint=0
+    evol_case=random.randint(N_observ//10,N_observ//5)
+    while coint < evol_case:
+        ch_obs=random.randint(0,len(SC[i])-1)
+        change=change_CN[random.randint(0,3)]
+        x_mean=(SC[i][ch_obs][2]*change)
+        SC[i][ch_obs][2]=int(random.normalvariate(x_mean, 2))
+        BSC[i][ch_obs][3]=SC[i][ch_obs][2]
+        a=random.random()
+        if a<=0.5:
+            if change==2:
+                BSC[i][ch_obs][2]=int(random.normalvariate(BSC[i][ch_obs][3]/2, 2))
+            elif change==1.5:
+                BSC[i][ch_obs][2]=int(random.normalvariate(BSC[i][ch_obs][3]/3, 2))
+            elif change==0:
+                BSC[i][ch_obs][2]=0
+            elif change==0.5:
+                BSC[i][ch_obs][2]=0
+                loh_stat=random.random()
+                if loh_stat<=0.1:
+                    SC[i][ch_obs][3]=int(random.normalvariate(SC[i][ch_obs][3]*2, 2))
+                    BSC[i][ch_obs][3]=SC[i][ch_obs][2]
+        coint+=1
+    VSC.append(BSC[i])
+    for j in range(len(VSC[i])):
+        mut=random.random()
+        if mut <=0.5:
+            if VSC[i][j][2]==0:
+                VSC[i][j][2]=int(random.normalvariate(VSC[i][j][3]/2,2))
+            else:
+                VSC[i][j][2]=int(random.normalvariate(VSC[i][j][2]/2,2))
+        else:
+            VSC[i][j][1]=int(random.normalvariate(VSC[i][j][1],N_observ/10))
+            VSC[i][j][2]=int(random.normalvariate(VSC[i][j][2]/2,2))
+
+        if VSC[i][j][2]<=0:
+            VSC[i][j][2]=1
+
+        if VSC[i][j][1] not in snv_dict.keys():
+            snv_dict[VSC[i][j][1]]=[[VSC[i][j][2]], VSC[i][j][3]]
+        else:
+            snv_dict[VSC[i][j][1]][0].append(VSC[i][j][2])
+
+print(len(snv_dict))
+
+CNA_data=copy.deepcopy(norm_data)
+for i in range(len(CNA_data)):
+    for j in SC:
+        CNA_data[i][2]+=j[i][2]
+    CNA_data[i][2]/=(len(SC)+1)
+    CNA_data[i][2]=int(CNA_data[i][2])
+
+BAF_data=copy.deepcopy(norm_data)
+for i in range(len(BAF_data)):
+    for j in BSC:
+        BAF_data[i][2]+=j[i][2]
+        BAF_data[i][3]+=j[i][3]
+    BAF_data[i][2]/=(len(BSC)+1)
+    BAF_data[i][2]=int(BAF_data[i][2])
+    BAF_data[i][3]/=(len(BSC)+1)
+    BAF_data[i][3]=int(BAF_data[i][3])
+
+# Output
+CNA=open('CNA_Emulate.txt', 'w')
+for i in CNA_data:
+    CNA.write(str(int(i[0]))+'\t'+str(int(i[1]))+'\t'+str(int(i[2]))+'\t'+str(int(i[3]))+'\n')
+CNA.close()
+
+BAF=open('BAF_Emulate.txt', 'w')
+for i in BAF_data:
+    BAF.write(str(int(i[0]))+'\t'+str(int(i[1]))+'\t'+str(int(i[2]))+'\t'+str(int(i[3]))+'\n')
+BAF.close()
+
+SNV=open('SNV_Emulate.txt', 'w')
+for i in snv_dict.keys():
+    snv_read=int(sum(snv_dict[i][0])/len(snv_dict[i][0]))
+    SNV.write('1'+'\t'+str(int(i))+'\t'+str(int(snv_read))+'\t'+str(int(snv_dict[i][1]))+'\n')
+SNV.close()
+
+
+print('-----end-----')
+#print(data)
+
+
+# In[ ]:
+
+
+