Mercurial > repos > yhoogstrate > varscan_mpileup2snp_from_bam

annotate test-data/generate_reads.py @ 91:94fb905e8d4f draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author yhoogstrate
date Thu, 05 Nov 2015 08:29:47 -0500
parents aaf98cdc5916
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
38
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
1 #!/usr/bin/env python
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
2
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
3
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
4 import random
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
5 import math
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
6
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
7
80
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
8 __version_info__ = ('1', '0', '0')
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
9 __version__ = '.'.join(__version_info__)
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
10
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
11
38
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
12 class Region:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
13 def __init__(self,start,stop,sequence):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
14 self.start = start
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
15 self.stop = stop
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
16 self.sequence = sequence.strip().replace("\n","").replace(" ","")
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
17 if(len(self.sequence) != self.getSpanningLength()):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
18 print "ERROR: sequence length: "+str(len(self.sequence))+", while spanning region is: "+str(self.getSpanningLength())
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
19 import sys
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
20 sys.exit()
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
21
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
22 def getSpanningLength(self):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
23 return abs(self.stop-self.start+1)
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
24
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
25 class ReadSynthesizer:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
26 def __init__(self,chromosome):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
27 self.regions = []
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
28 self.chromosome = chromosome
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
29
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
30 def addRegion(self,region):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
31 self.regions.append(region)
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
32
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
33 def produceReads(self,readDensity = 1,read_length = 50):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
34 """
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
35 Produces uniform reads by walking iteratively over self.regions
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
36 """
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
37
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
38 mRNA = self.getTotalmRNA()
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
39 spanning_length = self.getRegionSpanningLength()
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
40 n = spanning_length['total'] - read_length + 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
41
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
42 j = 0
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
43 k = 0
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
44
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
45 for i in range(n):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
46 # "alpha is playing the role of k and beta is playing the role of theta"
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
47 dd = max(0,int(round(random.lognormvariate(math.log(readDensity),0.5))))# Notice this is NOT a binomial distribution!!
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
48
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
49 for d in range(dd):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
50 sequence = mRNA[i:i+read_length]
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
51
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
52 if(random.randint(0,1) == 0):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
53 strand = 0
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
54 else:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
55 strand = 16
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
56 flag = strand + 0
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
57
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
58 print "read_"+str(j)+"."+str(i)+"."+str(d)+"\t"+str(flag)+"\t"+self.chromosome+"\t"+str(self.regions[j].start + k)+"\t60\t"+self.getMappingString(read_length,j,k)+"\t*\t0\t0\t"+str(sequence.upper())+"\t*"
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
59
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
60 spanning_length['iter'][j] -= 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
61 if(k >= self.regions[j].getSpanningLength()-1):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
62 j += 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
63 k = 0
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
64 else:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
65 k += 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
66
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
67 def getMappingString(self,length,j,offset):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
68 m = 0
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
69
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
70 out = ""
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
71
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
72 for i in range(length):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
73 k = i + offset
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
74
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
75 if(k >= self.regions[j].getSpanningLength()):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
76 j += 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
77
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
78 out += str(m)+"M"
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
79 out += (str(self.regions[j].start - self.regions[j-1].stop-1))+"N"
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
80 m = 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
81
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
82 offset = -k
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
83 else:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
84 m += 1
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
85
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
86 out += str(m) + "M"
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
87
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
88
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
89 return out
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
90
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
91 def getRegionSpanningLength(self):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
92 length = {'total':0,'iter':[]}
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
93 for r in self.regions:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
94 l = r.getSpanningLength()
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
95 length['iter'].append(l)
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
96 length['total'] += l
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
97 return length
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
98
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
99 def getTotalmRNA(self):
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
100 mRNA = ""
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
101 for r in self.regions:
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
102 mRNA += r.sequence
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
103 return mRNA
48c78adade03 Deleted selected files
yhoogstrate
parents:
diff changeset
104
80
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
105
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
106
5272ffe31456 Uploaded
yhoogstrate
parents: 38
diff changeset
107 if __name__ == "__main__":
88
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
108 # Real world example snp
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
109
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
110 #rs = ReadSynthesizer('chr6')
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
111 #rs.addRegion(Region(154360546,154360969,'ccaggactggtttctgtaagaaacagcaggagctgtggcagcggcgaaaggaagcggctgaggcgcttggaacccgaaaagtctcggtgctcctggctacctcgcacagcggtgcccgcccggccgtcagtaccatggacagcagcgctgcccccacgaacgccagcaattgcactgatgccttggcgtactcaagttgctccccagcacccagccccggttcctgggtcaacttgtcccacttagatggcGacctgtccgacccatgcggtccgaaccgcaccgacctgggcgggagagacagcctgtgccctccgaccggcagtccctccatgatcacggccatcacgatcatggccctctactccatcgtgtgcgtggtggggctcttcggaaacttcctggtcatgtatgtgattgtcag'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
112 #rs.addRegion(Region(154410961,154411313,'atacaccaagatgaagactgccaccaacatctacattttcaaccttgctctggcagatgccttagccaccagtaccctgcccttccagagtgtgaattacctaatgggaacatggccatttggaaccatcctttgcaagatagtgatctccatagattactataacatgttcaccagcatattcaccctctgcaccatgagtgttgatcgatacattgcagtctgccaccctgtcaaggccttagatttccgtactccccgaaatgccaaaattatcaatgtctgcaactggatcctctcttcagccattggtcttcctgtaatgttcatggctacaacaaaatacaggcaag'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
113 #rs.addRegion(Region(154412087,154412607,'gttccatagattgtacactaacattctctcatccaacctggtactgggaaaacctgctgaagatctgtgttttcatcttcgccttcattatgccagtgctcatcattaccgtgtgctatggactgatgatcttgcgcctcaagagtgtccgcatgctctctggctccaaagaaaaggacaggaatcttcgaaggatcaccaggatggtgctggtggtggtggctgtgttcatcgtctgctggactcccattcacatttacgtcatcattaaagccttggttacaatcccagaaactacgttccagactgtttcttggcacttctgcattgctctaggttacacaaacagctgcctcaacccagtcctttatgcatttctggatgaaaacttcaaacgatgcttcagagagttctgtatcccaacctcttccaacattgagcaacaaaactccactcgaattcgtcagaacactagagaccacccctccacggccaatacagtggatagaactaatcatcag'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
114 #rs.addRegion(Region(154428600,154428787,'gtggaattgaacctggactgtcactgtgaaaatgcaaagccttggccactgagctacaatgcagggcagtctccatttcccttcccaggaagagtctagagcattaattttgagtttgcaaaggcttgtaactatttcatatgatttttagagctgactatgacatgaaccctaaaattcctgttccc'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
115 #rs.produceReads(3,50)
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
116
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
117 # Artificial SNP
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
118 rs = ReadSynthesizer('chr1')
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
119 rs.addRegion(Region( 0+1, 59+1,'aaataggtcccaaacgttacgca'+'G'+'tctatgcctgacaaagttgcgaccacttcctctgcc'))#c -> G
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
120 rs.addRegion(Region( 60+1,119+1,'ttgtgtgacacgccggagatagg'+'A'+'catcagcaagtacgttaagtacactgaacgaactgg'))#g -> A
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
121 rs.addRegion(Region(120+1,179+1,'aggtttctacatcgtgcgtgatggc'+'C'+'ctaggagaagtgggtgtatctgcacagcataagt'))#t -> C
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
122 rs.addRegion(Region(180+1,239+1,'tataagacggaagtaaagcgtcttc'+'G'+'ccgttcagcaccccacgctcatagtcaatgctgg'))#a -> G
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
123 #rs.addRegion(Region(240+1,299+1,'ttcagcatagtcaagcgccggtggcctccaaaaagacgcactgagtagcttagctacttt'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
124 #rs.addRegion(Region(300+1,359+1,'gctccgcttgcggaagcactaagaggagattgaatttccaaatcccccccgatacctgtg'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
125 #rs.addRegion(Region(360+1,419+1,'cggtcgctacgtaagtgcgaagttctgttagatacgctccccttagtatatgggcgttaa'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
126 #rs.addRegion(Region(420+1,479+1,'tcggaccgtcggtactcactgcattccaggtctcatatagttcgccctagaagcctggga'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
127 rs.addRegion(Region(480+1,539+1,'tgaacgttgaacta'+'GCC'+'ctgatgtaaaccccgcgtgccaattccaggcgtcatgggggca'))#tag -> gcc
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
128 #rs.addRegion(Region(540+1,599+1,'acccctcgcagcctccctcttgctgttggtgcctagtatttcatgatttcgagccgacat'))
aaf98cdc5916 planemo upload for repository https://bitbucket.org/EMCbioinf/galaxy-tool-shed-tools commit 1477c39d48b290394b7247b9c7b1e4a62a85f2de-dirty
yhoogstrate
parents: 80
diff changeset
129 rs.produceReads(2,35)