Bioinformatics_scripts/random_reads.py at master · Dangertrip/Bioinformatics_scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
with open('/data/dsun/ref/mouseigenome/mm10.fa') as f:
    lines = f.readlines()
chr=''
ans={}
pos=0
mark=False
s=''
length=0
l=[]
chr_order=[]
for line in lines:
    temp = line.strip()
    if temp[0]=='>':
        if chr!='':
            ans[chr]=s
            ll=len(s)
            length+=ll
            l.append(ll)
            s=''
        chr= temp[1:]
        chr_order.append(chr)
        continue
    s+=temp
ans[chr]=s
ll=len(s)
l.append(ll)
length+=ll
import random
#print(ans)
def getcg(read):
    num=0
    for i in range(1,len(read)):
        if ('G'==read[i] or 'g'==read[i]) and ('C'==read[i-1] or 'c'==read[i-1]):
            num+=1
    return num

num=[]
print(l)
zero=0
for i in range(100000000):
    pos = random.randint(0,length)
    chr=0
    while pos>l[chr]:
        pos-=l[chr]
        chr+=1

    read=ans[chr_order[chr]][pos:min(pos+76,l[chr])]
    g=getcg(read)
    num.append(g)
    if g==0: zero+=1
print(zero/float(100000000))
print((100000000-zero)/float(100000000))

#import seaborn as sns
#from matplotlib import pyplot as plt
#import numpy as np
#n=np.array_split(num,10)
#for nn in n:
#    sns.distplot(nn)
#plt.savefig('125.pdf')


#