Perl随机产生多肽序列

yinyuemi

我的理解：
lz希望
1. 得到1000个的字符串，每个字符串的长度为7，组成这个字符串的单位（字符）是他由指定的20个字母
最终的数据类似如下
YGVDGQD
IVSKHMK
NYIWWAQ
WSHNLQW
SCPEWQT
...
2. 如果把上面的一系列字符串看成1000行 x 7列字母构成的话，每列中给字母出现的频率要符合下面的要求（即，每列的数值加起来为1）。
A        0.09497947        0.070426046        0.078231678        0.069402819        0.073001587        0.067751808        0.065225551
C        0.003987824        0.004239198        0.006126589        0.00579803        0.006359058        0.007059202        0.007186454
D        0.031654608        0.028931216        0.021204204        0.029173854        0.025034137        0.021646325        0.024539473
E        0.030191737        0.025704815        0.0171375               0.023008412        0.019931428        0.017206993        0.018367772
F        0.021186211        0.018803895        0.024108434        0.023349227        0.026832478        0.025921508        0.027780867
G        0.047137968        0.028274228        0.02906955        0.027688427        0.030174135        0.027163123        0.026097522
H        0.061280362        0.053520624        0.047650885        0.057602322        0.048347118        0.046961954        0.052401045
I        0.030525902        0.030684054        0.029671778        0.029735665        0.030828516        0.033772121        0.029383376
K        0.030053403        0.021923645        0.020619838        0.025112887        0.020833923        0.023215196        0.022188839
L        0.090544049        0.103690465        0.112389986        0.101604113        0.111333642        0.111322299        0.115992927
M        0.034562358        0.030696831        0.029966439        0.026875373        0.027945538        0.03101157        0.02620991
N        0.04517639        0.040286853        0.031861131        0.040137567        0.033102225        0.035602927        0.032904046
P        0.001273168        0.171672655        0.169522808        0.154487958        0.159433689        0.155568683        0.161045977
Q        0.064491638        0.052384357        0.043509473        0.051301545        0.043542068        0.040985172        0.047284907
R        0.03516902        0.033794677        0.046343819        0.04827906        0.055693288        0.058490476        0.059932226
S        0.136940411        0.097910949        0.106347363        0.10304508        0.104407819        0.105721405        0.103542222
T        0.133603708        0.103408973        0.098715659        0.093232499        0.088656527        0.099685822        0.083992158
V        0.05101914        0.041507999        0.040548266        0.040028439        0.043047795        0.039284095        0.041277095
W        0.022307746        0.013552552        0.018361253        0.015865506        0.019287738        0.018945619        0.018187325
Y        0.033914888        0.028585968        0.028613348        0.034271218        0.032207292        0.032683703        0.03646031

Monox

黎溪 发表于 2014-03-22 14:28 大家觉得我先产生10^9数量级的序列，再从里面随机抽取1000可行吗

从逻辑上是可行，不过既浪费空间又浪费时间，所以从实现上不一定可行（受到计算机资源的限制）。另一种方案是随机产生７个从０到１的数，这７个数分别按你指定的氨基酸７个位置的分布落入７个位置正好对应７个氨基酸这样就够成了一条７肽，重复１０００次就得到了１０００条序列，这个方案占用的空间可以不考虑（如果你是产生一条就输出一条到输出文件的话），时间也只与你要产生的序列的条数的线性相关（也就是O(n)），因此从实现上也是可行的。
代码还是希望你自己来写，当然就算我不帖代码，也多半会有其他人贴代码就是了，所以你也可以选择等待别人帖代码。不过，如果你不自己写，也不真正去理解别人给你的代码，那么别人给你的代码可能是错误的你也不知道。

黎溪

回复 11# yinyuemi
) 你的理解是正确的！

黎溪

回复 12# Monox
好的，谢谢你的建议！

   

Monox

我一直想学好Haskell，但是一直没有真正去用Haskell去完成什么实质性的任务，所以我就把你这个需求当成我自己的练习题了，希望以后尽量多用Haskell去解决实际问题。不过，你想要的代码是Perl的，我的代码可能对你没什么用。另一方面，希望有Haskell方面的专家能给我的代码提改进意见。好，下面是代码

1. $ cat amino.txt
   
2. A        0.09497947        0.070426046        0.078231678        0.069402819        0.073001587        0.067751808        0.065225551
   
3. C        0.003987824        0.004239198        0.006126589        0.00579803        0.006359058        0.007059202        0.007186454
   
4. D        0.031654608        0.028931216        0.021204204        0.029173854        0.025034137        0.021646325        0.024539473
   
5. E        0.030191737        0.025704815        0.0171375        0.023008412        0.019931428        0.017206993        0.018367772
   
6. F        0.021186211        0.018803895        0.024108434        0.023349227        0.026832478        0.025921508        0.027780867
   
7. G        0.047137968        0.028274228        0.02906955        0.027688427        0.030174135        0.027163123        0.026097522
   
8. H        0.061280362        0.053520624        0.047650885        0.057602322        0.048347118        0.046961954        0.052401045
   
9. I        0.030525902        0.030684054        0.029671778        0.029735665        0.030828516        0.033772121        0.029383376
   
10. K        0.030053403        0.021923645        0.020619838        0.025112887        0.020833923        0.023215196        0.022188839
    
11. L        0.090544049        0.103690465        0.112389986        0.101604113        0.111333642        0.111322299        0.115992927
    
12. M        0.034562358        0.030696831        0.029966439        0.026875373        0.027945538        0.03101157        0.02620991
    
13. N        0.04517639        0.040286853        0.031861131        0.040137567        0.033102225        0.035602927        0.032904046
    
14. P        0.001273168        0.171672655        0.169522808        0.154487958        0.159433689        0.155568683        0.161045977
    
15. Q        0.064491638        0.052384357        0.043509473        0.051301545        0.043542068        0.040985172        0.047284907
    
16. R        0.03516902        0.033794677        0.046343819        0.04827906        0.055693288        0.058490476        0.059932226
    
17. S        0.136940411        0.097910949        0.106347363        0.10304508        0.104407819        0.105721405        0.103542222
    
18. T        0.133603708        0.103408973        0.098715659        0.093232499        0.088656527        0.099685822        0.083992158
    
19. V        0.05101914        0.041507999        0.040548266        0.040028439        0.043047795        0.039284095        0.041277095
    
20. W        0.022307746        0.013552552        0.018361253        0.015865506        0.019287738        0.018945619        0.018187325
    
21. Y        0.033914888        0.028585968        0.028613348        0.034271218        0.032207292        0.032683703        0.03646031
    
22. $ cat gen_amino.hs
    
23. module Main where
    
24. 
    
25. import System.Environment (getArgs)
    
26. import System.Random (getStdGen, randoms)
    
27. 
    
28. parseAminoAcidDistribution :: String -> [(String, [Double])]
    
29. parseAminoAcidDistribution s = scanl1 f $ map fromRow rows
    
30.   where f (_, v) (k', v') = (k', zipWith (+) v v')
    
31.         fromRow row = (head row, map read (tail row))
    
32.         rows = map words $ lines s
    
33. 
    
34. genHeptapeptide :: [(String, [Double])] -> [Double] -> String
    
35. genHeptapeptide t ds = concatMap f [0..6]
    
36.   where f i = fst . head $ dropWhile (g i) t
    
37.         g i (_, v) | (ds !! i) >= (v !! i) = True
    
38.                    | otherwise = False
    
39. 
    
40. genHeptapeptides :: [(String, [Double])] -> [Double] -> [String]
    
41. genHeptapeptides t ds = genHeptapeptide t (take 7 ds) : genHeptapeptides t (drop 7 ds)
    
42. 
    
43. main = do
    
44.     args <- getArgs
    
45.     contents <- readFile "amino.txt"
    
46.     stdGen <- getStdGen
    
47.     let aminoAcidDist = parseAminoAcidDistribution contents
    
48.         randomNumbers = randoms stdGen
    
49.         n = read $ head args
    
50.         outFile = args !! 1
    
51.     writeFile outFile $ unlines $ take n $ genHeptapeptides aminoAcidDist randomNumbers
    
52. $ ghc gen_amino
    
53. [1 of 1] Compiling Main             ( gen_amino.hs, gen_amino.o )
    
54. Linking gen_amino ...
    
55. $ ./gen_amino 20 seq.txt
    
56. $ cat seq.txt
    
57. VATKDTV
    
58. ETDHDPP
    
59. TGPCSLS
    
60. SYLPDTH
    
61. TLKLSCN
    
62. HSSRPPP
    
63. KLLTSVL
    
64. YATKTPN
    
65. VPPLMFS
    
66. GPDMSGV
    
67. SPHLNWL
    
68. GKSPHVE
    
69. APHQKHS
    
70. QPADTLT
    
71. MAPTHVG
    
72. TSPNLSS
    
73. DVPVSLP
    
74. HPKALAL
    
75. THNLVPG
    
76. VGSLPSL
    

复制代码

从上面可以看到我写的代码是从命令行指定输出文件和输出序列的条数，我显示了输出２０条７肽的一次（因为结果是随机的，所以这个限定词很重要）运行结果。

pitonas

@Monox
@yinyuemi
师傅, 实在不好意思，比较茫然, 问题没看懂
哪位大神帮我理解这个问题?
能不能指导我一下呢？

Monox

pitonas 发表于 2014-03-25 15:31
@Monox
@yinyuemi
师傅, 实在不好意思，比较茫然, 问题没看懂

题外话，有人叫师傅还是蛮高兴的。
楼主的需求是根据一个七肽的分布矩阵随机生成１０００条七肽。
也就是说输入是一个七肽分布矩阵，输出是１０００条七肽。
这个分布矩阵有２０行７列。
所谓的七肽就是一个长度为７的字符串（准确来说是长度为７的氨基酸序列）。一般构成多肽（当然包括七肽）的是２０种氨基酸。输入的分布矩阵就是统计了７肽的７个位置的每种氨基酸的分布比例，然后产生的随机序列也得符合这个比例。
从另一个角度来说的话，也就是说假如我们产生的不是１０００条七肽，而是很多条（远大于１０００）七肽，理论上统计产生的所有七肽的各７个位置出现的每种氨基酸的次数会得到一个矩阵，这个矩阵的每个数除以总条数就应该得到一个分布矩阵，而这个分布矩阵应该和输入的矩阵完全一样。但是实际上这种统计学的问题是肯定不会完全一样的，肯定会有一些偏差，因为没办法产生无数条序列来进行统计只能取样统计嘛。
我担心没说明白，就再说一遍，具体那个分布矩阵是怎么来的。其实一般那个分布矩阵是根据已有的七肽序列库先进行了统计而得到的，比如说其中有两条七肽，一条为“VATKDTV”而另一条为“ETDHDPP”，那么首先在Ｖ所对应的行的第一列加１（当然初始值是０），在Ａ所对应的行的第二列加１，Ｔ所对应的行的第三列加１……，在Ｅ所对应的行的第一列加１，Ｔ所对应的行的第二列加１，Ｄ所对应的行的第三列加１……，对每条七肽都这样处理以后得到的一个２０行７列的矩阵里的每个数除以总条数就是作为楼主问题的输入的分布矩阵。而这个分布矩阵的一个特点是每列的和近似等于１（因为有浮点误差，所以一般不可能正好相等了，楼主给的输入矩阵的每列的和就稍高于１一点点）。

ba_du_co

表示没看懂

yinyuemi

回复 16# pitonas


   我想 Monox大神已经说得很清楚了~
   
   题外话，越来越觉得，想在perl版混的话，不仅要懂perl，还得掌握一些的生物学知识，这里快被生物信息学童鞋们的题目侵占了。
  既来之，则安之。
  所以，建议童鞋们发帖之前，最好能把生物学学语言翻译成计算机语言才能让你的问题得到更好地解决，让外行听懂你的意图，这也是一种能力，童鞋们要加油啊~
   

Monox

回复 18# ba_du_co
既然还有人说没看懂，我就再来讲个故事吧，用来说明这个问题。因为有人叫师傅所以才这么耐心的，要是还有人说看不懂，就只能说我当不了师傅，就等当得了师傅的人再来解释

有一个生物信息学的师傅给徒弟安排了一个任务，任务很简单，师傅给了徒弟两条3肽序列，一条是AEY，另一条是DEC，师傅让徒弟统计组成这两条3肽的二十种氨基酸是如何分布的。这个徒弟还说这么简单的问题，他首先在一个EXCEL表格里前二十行的第一列写上二十种氨基酸的单字母缩写形式，并而他把前6行按A，B，C，D，E,Y这个形式写，还有14行其他的字母，他首先看一下AEY，很简单，他往EXCEL里填入一些数字，结果成了下面的形式(列出前七行，后面的行省略)

1. 
   
2. A 1 0 0
   
3. B 0 0 0
   
4. C 0 0 0
   
5. D 0 0 0
   
6. E 0 1 0
   
7. Y 0 0 1
   
8. P 0 0 0
   

复制代码

他再看了一眼下条序列，DEC，结果他把EXCEL改成了下面的这个样子（同样只列出前七行）：

1. 
   
2. A 1 0 0
   
3. B 0 0 0
   
4. C 0 0 1
   
5. D 1 0 0
   
6. E 0 2 0
   
7. Y 0 0 1
   
8. P 0 0 0
   

复制代码

他把这个结果给了他师傅。他师傅一看，“嗯，很好，不过这个分布做一下归一化吧，这个结果和其他的数据整合不起来”。于是徒弟知道师傅的意思是把每个数除以２（序列的总条数为２），徒弟把修改好的结果给了师傅，结果的前七行是这样子的：

1. 
   
2. A 0.5 0 0
   
3. B 0    0 0
   
4. C 0    0 0.5
   
5. D 0.5 0 0
   
6. E 0    1 0
   
7. Y 0   0 0.5
   
8. P 0   0 0
   

复制代码

师傅提醒徒弟，每一列的和应该是1，这样可以用来检验他的结果有没有计算差误，徒弟按师傅的提醒统计了一下，果然每一列的和都是1。师傅接着说，“你的第一个问题完成得很好，我现在要你根据你得到的这个分布表随机生成3条3肽”。徒弟有点犯难，他不知道如何随机生成3肽，还好，他找到了一个程序可以达到这样的效果。这个程序首先产生了一个字符，这个字符为A，这个程序产生另一个字符是E，这个程序产生第三个字符C，于是第一条3肽出现了。同样的方式，第二条3肽也出现了，接着是第三条。这个徒弟虽然不会写程序，但是他知道产生的3肽有如下一些特点，就是第二个氨基酸为E的概率为1（分布表里的数），而为其他氨基酸的概率为０，也就是说产生的所有的序列的第二个位置都是E，而同样的道理，第一个位置只会出现A和D，而它们出现的概率相同，而第三个位置只能出现C和Y，出现的概率同样都是50%。徒弟按这个特点检查了结果后把结果给师傅了，还把这些特点也跟师傅说了。师傅说，“很好，你用这3条3肽再按问题1统计一下，看有什么规律。再生成10条3肽再统计一下看看，再生成1000条看看。”徒弟按师傅的要求做了，他发现用生成的3肽统计的结果本来和原来的EXCEL的结果差别还是蛮大的（除了为0和为1的格子没有变化外），但是随着生成的序列的条数不断增加，徒弟发现统计结果和原来的EXCEL的结果越来越接近了。师傅很高兴徒弟告诉了他这个结论，于是师傅说，“很好，我这里有500条七肽，我按给你的问题一的方式已经统计过氨基酸的分布了，你就直接根据那个分布随机生成1000条七肽吧，相信你在处理完3肽的问题后可以做得来的”。而神奇的是师傅给的那个分布就是楼主给的那个矩阵，徒弟知道结果应该满足的特点，但是他不知道如何修改那个程序才能让程序产生七肽而不是三肽，因为徒弟完全不会计算机，那么各位有可以帮他改写那个程序的吗？喔，不好，徒弟电脑种病毒那个程序被删了，所以各位好心的网友只能帮徒弟重写而不是改原来的程序了。于是我就写了上面的Haskell的程序给那个徒弟，那个徒弟很高兴，他的师傅也很高兴，不过楼主不怎么高兴。因为楼主的问题虽然和那个师傅给徒弟的问题一样，但是楼主只要Perl的代码，他不想发那么多时间去下载那么庞大的Haskell的编译器来用我写的Haskell程序，而且楼主说Perl是世界上最好的语言，他只要Perl的答案。所以就请各位Perl高手帮帮楼主解决这个问题吧。

总之，楼主给的矩阵里的数字是要产生的七肽的七个位置不同氨基酸出现的概率。