Linux脚本,以获得所有可能的7个字母组合,以产生肽在pymol



我希望生成一个文件夹,包含7(长度)特定氨基酸的每个肽的pdb文件。我想首先制作一个简单的linux脚本来生成一个文件,其中包含所有7个字母的组合,如:

AAAAAAA
AAAAAAB
AAAAABA
AAAABAA
AAABAAA
AABAAAA
ABAAAAA
BAAAAAA
AAAAABB
AAAABAB
...

我认为这个脚本可以工作,但我不确定:

for c1 in {A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
do
    for c2 in {A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
    do
        for c3 in {A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
        do
            for c4 in {A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
            do
                for c5 in {A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
                do
                    printf "%sn" "$c1$c2$c3$c4$c5"
                done
            done
        done
    done
done

然后使用其他简单的脚本,其中最后一个文件的每一行都使用pymol生成肽,并使用以下命令:

for aa in "row1": cmd._alt(string.lower(aa))
save row1.pdb, all

我是linux脚本新手。请问有人能帮我吗?谢谢

我看了看(ab?)的想法,使用大括号展开:

p='{A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}'
eval echo $p$p$p$p$p$p$p

在7个$p的简单步骤中使用这种直接的方法对bash来说太过分了。没有明显的原因,它吞噬了所有的内存(随时间的测量显示没有其他内存值增加得如此之快)。对于多达4个$p,该命令非常快速和简单,只有两行:

p='{A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}'
eval echo $p$p$p$p
然而,内存使用增长非常快。在6个$p重复的深度,该过程消耗超过7.80 gb的内存。eval部分还有助于增加执行时间和内存使用。

需要另一种方法。所以,我试着让每一步的扩张都独立进行,利用乔纳森·莱弗勒使用的概念。对于输入中的每一行,编写19行,每一行在输出中添加一个字母。我发现任何eval都是一个重要的内存消耗(这里没有显示)。

Bash

更简单的bash过滤器是:

bashfilter(){
    while read -r line; do
        printf '%sn' ${line}{A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
    done </dev/stdin
}

可用于以下几个级别的处理:

echo | bashfilter | bashfilter | bashfilter

它只需要按照每行需要的字母重复过滤器步骤。

使用这个更简单的方法:内存不再是一个问题。然而,速度变差了。

Leffler SED

只是为了比较,用它作为一个标尺,我实现了莱弗勒的想法:

# Building Leffler solution:
    leftext="$(<<<"${list}" sed -e 's/,/n/g')"                 # list into a column.
    leftext="$(<<<"${leftext}" sed -e 's%.%s/$/&/p;s/&$//%')"   # each line ==> s/$/?/p;s/?$//
    # echo -e "This is the leffilter n$leftext"
leffilter(){ sed -ne "$leftext"; }    # Define a function for easy use.

And是可以递归地使用的字符过滤器,可以根据需要获得每行任意多的字母:

echo | leffilter | leffilter | leffilter

Leffler解决方案插入一个字母,删除一个字母。


SED

不需要擦去一个字母,工作量就可以减少。我们可以将原始模式空间存储在"保持空间"中。

然后,将第一行复制到保持空间(h),并继续恢复(g)并插入一个字母。

# Building a sed solution:
    sedtext="$(<<<"${list}" sed -e 's/,/n/g')"    # list into a column.
    sedtext="$(<<<"${sedtext}" sed -e 's%[A-Z]%g;s/$/&/p;%g')"  # s/$/?/p
    sedtext="$(<<<"${sedtext}" sed -e '1 s/g/h/' )"             # 1st is h
sedfilter(){ sed -ne "$sedtext"; }    # Define a function for easy use.  

这样做可以提高速度,大约降低1/3(33%)。或者快1.47倍。


AWK

最后,我提出了一个AWK解决方案。我之前写过,但它是最快的。所以我把它作为最后一个选项。最好的,直到有人提出更好的:-)

# An AWK based solution:
awkfilter(){ awk 'BEGIN { split( "'"$list"'",l,",");}
                        { for (i in l) print $0 l[i] }'
}

对,只有两行。它的速度是Leffler溶液的一半或两倍。

使用的完整测试脚本如下。它重新调用自身以启用外部时间的使用。确保它是一个使用bash的可执行文件。

#!/bin/bash
TIMEFORMAT='%3lR %3lU %3lS'
list="A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y"
# A pure bash based solution:
bashfilter(){
    while read -r line; do
        printf '%sn' ${line}{A,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
    done </dev/stdin
}
# Building Leffler solution:
    leftext="$(<<<"${list}" sed -e 's/,/n/g')"                 # list into a column.
    leftext="$(<<<"${leftext}" sed -e 's%.%s/$/&/p;s/&$//%')"   # each line ==> s/$/?/p;s/?$//
    # echo -e "This is the lef filter n$leftext"
leffilter(){ sed -ne "$leftext"; }    # Define a function for easy use.
# Building a sed solution:
    sedtext="$(<<<"${list}" sed -e 's/,/n/g')"                 # list into a column.
    sedtext="$(<<<"${sedtext}" sed -e 's%[A-Z]%g;s/$/&/p;%g')"  # each letter ==> s/$/?/p
    sedtext="$(<<<"${sedtext}" sed -e '1 s/g/h/' )"             # First command is 'h'.
    # echo -e "This is the sed filter n$sedtext"
sedfilter(){ sed -ne "$sedtext"; }    # Define a function for easy use.
# An AWK based solution:
awkfilter(){ awk 'BEGIN { split( "'"$list"'",l,",");}
                        { for (i in l) print $0 l[i] }'
}
# Execute command filter
docommand(){
    local a count="$1" filter="$2" peptfile="$3"
    for (( i=0; i<count; i++ )); do
        case $filter in
            firsttry) a+=("{$list}"); ;;
            *)        a+=("| $filter"); ;;
        esac
    done
    [[ $filter == firsttry ]] && a+=('| sed '"'"'s/ /n/'"'" )
    [[ -n $peptfile ]] && peptfile="$peptfile.$count"
    eval 'echo '"$(printf '%s' "${a[@]}")" > "${peptfile:-/dev/null}";
}
callcmd(){
    tf='wall:%e s:%S u:%U (%Xtext+%Ddata %F %p %t %Kmem %Mmax)'
    printf '%-12.12s' "$1" >&2
    /usr/bin/time -f "$tf" "$0" "$repeats" "$1" "$2"
}
nofile=1
if (( $#>=2 )); then
    docommand "$1" "$2" "$3"; exit 0
else
    for (( i=1; i<=6; i++)); do
        repeats=$i; echo "repeats done = $repeats"
        if ((nofile)); then
            callcmd firsttry
            callcmd bashfilter
            callcmd leffilter
            callcmd sedfilter
            callcmd awkfilter
        else
            callcmd firsttry   peptidesF
            callcmd bashfilter peptidesB
            callcmd leffilter  peptidesL
            callcmd sedfilter  peptidesS
            callcmd awkfilter  peptidesA
        fi
    done
fi
<标题> 结果

使用外部程序/usr/bin/time(而不是bash内置时间)来测量所使用的内存。这在这个问题中很重要。

: tf = '墙:% e s: % s u: % u (% Xtext + % Ddata % F % p % t % Kmem % Mmax)"

使用上面的脚本很容易找到7个循环和真实文件输出的结果,但是我觉得填充大约21 gb的磁盘空间太多了。

6次循环前的结果如下:

   repeats done = 1
firsttry    wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)
bashfilter  wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1552max)
leffilter   wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)
sedfilter   wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)
awkfilter   wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1560max)

:

   repeats done = 2
firsttry    wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)
bashfilter  wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1552max)
leffilter   wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1560max)
sedfilter   wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)
awkfilter   wall:0.01 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1560max)

:

   repeats done = 3
firsttry    wall:0.02 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1796max)
bashfilter  wall:0.07 s:0.00 u:0.05 (0text+0data 0 0 0 0mem 1552max)
leffilter   wall:0.02 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)
sedfilter   wall:0.02 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1560max)
awkfilter   wall:0.02 s:0.00 u:0.00 (0text+0data 0 0 0 0mem 1556max)

:

   repeats done = 4
firsttry    wall:0.28 s:0.01 u:0.26 (0text+0data 0 0 0 0mem 25268max)
bashfilter  wall:0.96 s:0.03 u:0.94 (0text+0data 0 0 0 0mem 1552max)
leffilter   wall:0.13 s:0.00 u:0.12 (0text+0data 0 0 0 0mem 1560max)
sedfilter   wall:0.10 s:0.00 u:0.08 (0text+0data 0 0 0 0mem 1560max)
awkfilter   wall:0.09 s:0.00 u:0.07 (0text+0data 0 0 0 0mem 1560max)

:

   repeats done = 5
firsttry    wall:4.98 s:0.36 u:4.76 (0text+0data 0 0 0 0mem 465100max)
bashfilter  wall:20.19 s:0.81 u:20.18 (0text+0data 0 0 0 0mem 1552max)
leffilter   wall:2.43 s:0.00 u:2.50 (0text+0data 0 0 0 0mem 1556max)
sedfilter   wall:1.83 s:0.01 u:1.87 (0text+0data 0 0 0 0mem 1556max)
awkfilter   wall:1.49 s:0.00 u:1.54 (0text+0data 0 0 0 0mem 1560max)

:

   repeats done = 6
firsttry    wall:893.06 s:30.04 u:105.22 (0text+0data 402288 0 0 0mem 7802372m)
bashfilter  wall:365.13 s:14.95 u:368.09 (0text+0data 0 0 0 0mem 1548max)
leffilter   wall:51.90 s:0.09 u:53.91 (0text+0data 6 0 0 0mem 1560max)
sedfilter   wall:35.17 s:0.08 u:36.67 (0text+0data 0 0 0 0mem 1556max)
awkfilter   wall:25.60 s:0.06 u:26.77 (0text+0data 1 0 0 0mem 1556max)

这里有一个技巧可以"相当快"地得到答案。基本上,它从一个包含单个换行符和氨基酸字母列表的文件开始。它生成一个sed脚本(当然是使用sed),它连续地在一行的末尾添加一个氨基酸字母,打印它,删除它,然后移动到下一个字母。

peptides-A.sh

printf "%sn" A D E F G H I K L M N P Q R S T V W Y |
sed 's%.%s/$/&/p;s/&$//%' > peptides.sed
echo > peptides.0A      # Bootstrap the process
        sed -n -f peptides.sed peptides.0A > peptides.1A
        sed -n -f peptides.sed peptides.1A > peptides.2A
        sed -n -f peptides.sed peptides.2A > peptides.3A
timecmd sed -n -f peptides.sed peptides.3A > peptides.4A
timecmd sed -n -f peptides.sed peptides.4A > peptides.5A
timecmd sed -n -f peptides.sed peptides.5A > peptides.6A
timecmd sed -n -f peptides.sed peptides.6A > peptides.7A

你可以把'timecmd'看作是time的一个变体。它打印开始时间、命令,然后运行它,然后打印结束时间和经过的时间(仅限时钟时间)。

样本输出:

$ bash peptides-A.sh
2015-10-16 15:25:24
+ exec sed -n -f peptides.sed peptides.3A
2015-10-16 15:25:24 - elapsed: 00 00 00
2015-10-16 15:25:24
+ exec sed -n -f peptides.sed peptides.4A
2015-10-16 15:25:27 - elapsed: 00 00 03
2015-10-16 15:25:27
+ exec sed -n -f peptides.sed peptides.5A
2015-10-16 15:26:16 - elapsed: 00 00 49
2015-10-16 15:26:16
+ exec sed -n -f peptides.sed peptides.6A
2015-10-16 15:42:47 - elapsed: 00 16 31
$ ls -l peptides.?A; rm -f peptides-?A
-rw-r--r--  1 jleffler  staff           1 Oct 16 15:25 peptides.0A
-rw-r--r--  1 jleffler  staff          38 Oct 16 15:25 peptides.1A
-rw-r--r--  1 jleffler  staff        1083 Oct 16 15:25 peptides.2A
-rw-r--r--  1 jleffler  staff       27436 Oct 16 15:25 peptides.3A
-rw-r--r--  1 jleffler  staff      651605 Oct 16 15:25 peptides.4A
-rw-r--r--  1 jleffler  staff    14856594 Oct 16 15:25 peptides.5A
-rw-r--r--  1 jleffler  staff   329321167 Oct 16 15:26 peptides.6A
-rw-r--r--  1 jleffler  staff  7150973912 Oct 16 15:42 peptides.7A
$

我使用问题中的脚本创建peptides.5B(脚本在我的磁盘上被称为peptides-B.sh),并检查peptides.5Apeptides.5B是相同的。

测试环境:13" MacBook Pro, 2.7 GHz Intel Core i5, 8gb RAM, SSD存储。


编辑行首而不是行尾可使性能提高约20%。

代码:

printf "%sn" A D E F G H I K L M N P Q R S T V W Y |
sed 's%.%s/^/&/p;s/^&//%' > peptides.sed
echo > peptides.0A      # Bootstrap the process
        sed -n -f peptides.sed peptides.0A > peptides.1A
        sed -n -f peptides.sed peptides.1A > peptides.2A
        sed -n -f peptides.sed peptides.2A > peptides.3A
timecmd sed -n -f peptides.sed peptides.3A > peptides.4A
timecmd sed -n -f peptides.sed peptides.4A > peptides.5A
timecmd sed -n -f peptides.sed peptides.5A > peptides.6A
timecmd sed -n -f peptides.sed peptides.6A > peptides.7A

时间:

$ bash peptides-A.sh; ls -l peptides.?A; wc peptides.?A; rm -f peptides.?A
2015-10-16 16:05:48
+ exec sed -n -f peptides.sed peptides.3A
2015-10-16 16:05:48 - elapsed: 00 00 00
2015-10-16 16:05:48
+ exec sed -n -f peptides.sed peptides.4A
2015-10-16 16:05:50 - elapsed: 00 00 02
2015-10-16 16:05:50
+ exec sed -n -f peptides.sed peptides.5A
2015-10-16 16:06:28 - elapsed: 00 00 38
2015-10-16 16:06:28
+ exec sed -n -f peptides.sed peptides.6A
2015-10-16 16:18:51 - elapsed: 00 12 23
-rw-r--r--  1 jleffler  staff           1 Oct 16 16:05 peptides.0A
-rw-r--r--  1 jleffler  staff          38 Oct 16 16:05 peptides.1A
-rw-r--r--  1 jleffler  staff        1083 Oct 16 16:05 peptides.2A
-rw-r--r--  1 jleffler  staff       27436 Oct 16 16:05 peptides.3A
-rw-r--r--  1 jleffler  staff      651605 Oct 16 16:05 peptides.4A
-rw-r--r--  1 jleffler  staff    14856594 Oct 16 16:05 peptides.5A
-rw-r--r--  1 jleffler  staff   329321167 Oct 16 16:06 peptides.6A
-rw-r--r--  1 jleffler  staff  7150973912 Oct 16 16:18 peptides.7A
        1         0          1 peptides.0A
       19        19         38 peptides.1A
      361       361       1083 peptides.2A
     6859      6859      27436 peptides.3A
   130321    130321     651605 peptides.4A
  2476099   2476099   14856594 peptides.5A
 47045881  47045881  329321167 peptides.6A
893871739 893871739 7150973912 peptides.7A
943531280 943531279 7495831836 total
$

我启动了wc的输出,所以它是"正确列"(添加空格,换句话说)。当数字包含8位时,原来的数字开始不稳定

免责声明:虽然我很高兴已经找到了这个算法基于基数19的数字,它是难以忍受的慢(8秒3个字母的字符串,160秒4个字母的,都有19个氨基酸,运行在2.2 GHz核心i7没有实际节省输出)相比其他解决方案,乔纳森·勒夫勒暗示。不管怎样,我还是把它留在这里,以防有人和我一样觉得有趣。

这是一个可能的替代方案,最多有19个氨基酸(你在代码中引用的那些):

aminoarr=("A" "D" "E" "F" "G" "H" "I" "K" "L" "M" "N" "P" "Q" "R" "S" "T" "V" "W" "Y")
peplength=7
aminonum=19
N=0
while [ $N -le $(( ${aminonum}**${peplength} - 1 )) ]; do
  remain=$N
  #printf "%d " $N
  for k in $(seq $(( ${peplength}-1 )) -1 0 ) ; do
    digit=$(( ${remain} / (${aminonum}**${k}) ))
    printf "%s" ${aminoarr[$digit]}
    let remain=$(( ${remain} - ${digit}*(${aminonum}**${k}) ))
  done
  echo
  let N=${N}+1
done

首先我们定义了氨基酸序列(aminoarr),我们不能生成的肽的长度(peplength),以及我们想要从列表中选择的氨基酸数量(aminonum,不应大于19)。

然后我们从N循环到aminonum^peplength -1,基本上生成以19为基数的所有可能的数字,最多7位(如果我们坚持你的问题中的参数)。然后将每个数字以19为基数进行分解,并从阵列aminoarr中选择相应的氨基酸。请注意,在基数19中,每个数字都落在0到18之间,因此它们非常适合索引包含19个元素的aminoarr

如果你取消注释printf行,它会给你给定序列的数字,但这会使你的文件更大(正如@Jonathan Leffler非常正确地评论了输出大小)。

无论如何,下面是前20行的示例输出:
AAAAAAA
AAAAAAD
AAAAAAE
AAAAAAF
AAAAAAG
AAAAAAH
AAAAAAI
AAAAAAK
AAAAAAL
AAAAAAM
AAAAAAN
AAAAAAP
AAAAAAQ
AAAAAAR
AAAAAAS
AAAAAAT
AAAAAAV
AAAAAAW
AAAAAAY
AAAAADA

crunch在Kali发行版上可用

crunch 7 7 ADEFGHIKLMNPQRSTVWY

最新更新