Jun 18 2008

爱的方程

分类: programming — azalea @ 5:36 am 250次阅读

爱的方程： 17x^2-16|x|y+17y^2=225

先从风雷的技术天地看到：http://www.forwind.cn/2008/06/09/draw-heart/

又寻根溯源找到Matrix67：http://www.matrix67.com/blog/archives/85

两个都是很赞的blog，推荐一下！（汗，跑题了。。）

原理很简单，就是画出方程 17x^2-16|x|y+17y^2=225  的图像：

以下是代码（基本都是从风雷的技术天地那里抄来的，不过窃以为我的配色比较正常@@）

PLAIN TEXT
PYTHON:
  1. from pylab import *
  2.  
  3. def func1(x):
  4.     return ( 16*fabs(x) + sqrt( 15300 - (x**2) * 900 ) ) / 34
  5. def func2(x):
  6.     return ( 16*fabs(x) - sqrt( 15300 - (x**2) * 900 ) )/ 34
  7.  
  8. x = arange( -5.0,5.0,0.001 )
  9. plot(x,func1(x),c='r')
  10. plot(x,func2(x),c='r')
  11. show()

以上代码需要matplotlib，Ubuntu可以sudo apt-get install python-matplotlib

另外这里还有一个更高级的方程，不过比较丑陋。

相关文章


Jun 05 2008

使用python进行ClustalW格式和aligned FASTA格式转换

分类: programming, study — azalea @ 5:49 pm 311次阅读

ClustalW格式是多序列比对程序CLUSTAL W的默认输出格式，大概是这个样子：

Aligned FASTA格式是多序列比对程序MUSCLE的默认输出格式，与一般的FASTA格式很相似，只是根据比对结果在序列中填充“-”，使所有序列长度相同，例如：

Aligned FASTA格式更便于机器阅读，而ClustalW格式是比较人性化的 ，使用者可以直观地看到相同和不同的区域。

由于MUSCLE程序的profile alignment功能需要aligned FASTA格式的输入，因此我写了以下一段简单的python代码，把ClustalW格式转换成aligned FASTA格式：

PLAIN TEXT
PYTHON:
  1. #!/usr/bin/python
  2. #Usage: python clw2afa.py input.clw
  3. #Default output is STDOUT, you can redirect using '>'
  4.  
  5. import sys
  6. import re
  7.  
  8. pat=re.compile('(\S+)\s+(\S+)')
  9.  
  10. f=open(sys.argv[1])
  11. header=f.readline()
  12. if not header.startswith('CLUSTAL') and not header.startswith('MUSCLE'):
  13.     sys.exit("Wrong input file type.")
  14. ln=f.readline()
  15. while not ln.strip():
  16.     ln=f.readline()
  17. seqdict={}
  18. EOF=False
  19. while ln:
  20.     m=pat.search(ln)   
  21.     if m and not '*' in ln:
  22.         seqdict[m.group(1)]=m.group(2)
  23.         ln=f.readline()
  24.     elif not ln:
  25.         EOF=True
  26.         break
  27.     else:
  28.         break
  29. if not EOF:
  30.     while True:
  31.         ln=f.readline()
  32.         if not ln:
  33.             break
  34.         m=pat.search(ln)
  35.         if m and not '*' in ln:
  36.             seqdict[m.group(1)]+='\n'+m.group(2)
  37. f.close()
  38. for k,v in seqdict.items():
  39.     print ">%s\n%s"%(k,v)

输入文件

写得有点混乱，希望大家如果使用的话把BUG报告一下^^

相关文章


May 15 2008

Python下的lowess拟合

分类: programming — azalea @ 2:27 pm 254次阅读

前段时间用R的lowess function用得很爽，于是想用python也lowess一下。找到了这个东东：

http://mail.python.org/pipermail/python-list/2004-September/281365.html

发现竟然是我崇拜的Istvan写的，是把R里lowess function的C代码包装成python代码。

下载后解压，安装：

Linux下：
sudo apt-get install swig #因为这个module依赖于swig (Generate scripting interfaces to C/C++ code)
python setup.py build
python setup.py install

Windows下：
把win32/2.3 子目录下的2个文件copy到python/site_packages目录

用法：
import lowess
ys = lowess.fit(x=X, y=Y, F=F, NSTEPS=NSTEPS, DELTA=DELTA)
具体请参考解压后文件夹中的 lowess_readme.txt

举例：

PLAIN TEXT
PYTHON:
  1. from pylab import *
  2. import lowess
  3. # 需要matplotlib module, 详见http://matplotlib.sourceforge.net/
  4. length=[176, 218, 113, 102, 250, 309, 135, 150, 200, 263, 135, 114, 123, 118, 100, 109, 100, 94, 184, 303, 117, 80, 110, 108, 123, 85, 85, 138, 163, 184, 196, 225, 172, 199, 206, 118, 108, 98, 123, 99, 125, 124, 133, 95, 102, 93, 102, 85, 89, 89, 237, 236, 331, 173, 85, 257, 331, 111, 107, 120, 117, 116, 123, 111, 140, 221, 182, 111, 337, 191, 178, 177, 191, 205, 186, 135, 151, 195, 400, 249, 264, 207, 113, 153, 102, 208, 184, 107, 109, 187, 187, 333, 138, 160, 109, 377, 337, 213, 221, 206, 190, 215, 211, 411, 181, 226, 379, 167, 217, 154, 106, 73, 193, 295, 118, 117, 121, 168, 138, 124, 689, 163, 105, 122, 113, 143, 120, 139, 170, 198, 112, 100, 95, 100, 100, 93, 126, 92, 142, 103, 126, 132, 212, 221, 86, 124, 157, 174, 189, 133, 160, 111, 100, 106, 123, 104, 183, 107, 118, 103, 132, 226, 212, 199, 185, 359, 230, 195, 205, 221, 105, 115, 116, 221, 221, 377, 271, 272, 290, 126, 221, 132, 217, 98]
  5. n, bins, patches = hist(length, 100, normed=0)
  6. # Draw a histogram with 100 bins
  7. # n is a list of frequency of each bin; bin is a list of the leftmost position of bins
  8. # hist()具体参数详见http://matplotlib.sourceforge.net/matplotlib.pyplot.html#-hist
  9. ys = lowess.fit(x=list(bins), y=list(n), F=0.2, NSTEPS=2, DELTA=0.2)
  10. plot(list(bins),ys,c='r',linewidth=2)
  11. # plot()具体参数详见http://matplotlib.sourceforge.net/matplotlib.pyplot.html#-plot
  12. show()

图：
lowess

此外，Biopython的Statistics包也有lowess模块，不会用-,-

相关文章


May 15 2008

Python计算中位数和众数

分类: programming — azalea @ 11:38 am 356次阅读

一些代码：

PLAIN TEXT
PYTHON:
  1. def median(numbers):
  2. '''Return the median of the list of numbers.'''
  3. # Sort the list of numbers and take the middle element.
  4. n = len(numbers)
  5. copy = numbers[:] # So that "numbers" keeps its original order
  6. copy.sort()
  7. if n & 1:         # There is an odd number of elements
  8. return copy[n / 2]
  9. else:
  10. return (copy[n / 2 - 1] + copy[n / 2]) / 2
  11.  
  12. def mode(numbers):
  13. '''Return the mode of the list of numbers.'''
  14. #Find the value that occurs the most frequently in a data set
  15. freq={}
  16. for i in range(len(numbers)):
  17. try:
  18. freq[numbers[i]] += 1
  19. except KeyError:
  20. freq[numbers[i]] = 1
  21. max = 0
  22. mode = None
  23. for k, v in freq.iteritems():
  24. if v> max:
  25. max = v
  26. mode = k
  27. return mode

同样的事用python stats module也可以做，对应的函数就是stats.median()和stats.mode()，但是我用stats.median()结果很奇怪，输入是一堆整数的list，结果median是小数。。希望高手告诉我这是为啥。。stats.mode返回一个list，第一个元素是众数的频度，第二个元素是所有众数的list（因为未必唯一），而上面的代码仅返回其中一个众数，当然简单改进一下就和stats.mode()完全一样啦

Reference:
http://mail.python.org/pipermail/python-list/2004-December/294990.html

相关文章


Apr 16 2008

python函数的问题呀。。

分类: programming — azalea @ 2:02 pm 282次阅读

没时间查暂时，记下来先

如果一个函数的目的是改变字典中一个key的值，我是传递这个值进去，还是传递整个字典进去？

如果我传递字典进去，是引用传递还是值传递？在函数里改变key的值， 如果没有返回值，那么原来的字典到底有没有改变？

如果不改变，那么怎么才能引用传递参数，或者说传递参数的地址给函数？

答案：是引用传递啦

相关文章


第 1 页，共 4 页1234>