0

我正在寻找一个包含 7 个(lentgh)特定氨基酸的每个肽的 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 "%s\n" "$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 脚本的新手。有人可以帮助我吗?谢谢

4

4 回答 4

3

我看了一下(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 Gigs 的内存。eval 部分还有助于增加执行时间和内存使用量。

需要一种替代方法。因此,我尝试利用 Jonathan Leffler 使用的概念,自行完成扩展的每一步。对于输入中的每一行,写 19 行,每行在输出中附加一个字母。我发现任何 eval 都是重要的内存消耗(此处未显示)。

重击

一个更简单的 bash 过滤器是:

bashfilter(){
    while read -r line; do
        printf '%s\n' ${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

它只需要重复尽可能多的过滤步骤,就像每行需要的字母一样。

使用这种更简单的方法:内存不再是问题。然而,速度变得更糟。

莱弗勒 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.

并且是可以递归使用的 leffilter 以获得每行所需的尽可能多的字母:

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 '%s\n' ${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 内置时间)来测量所使用的内存。在这个问题中很重要。

With: tf='wall:%es:%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)
于 2015-10-19T21:02:33.297 回答
1

这是一种产生“相当快”答案的技术。基本上,它从一个包含单个换行符和氨基酸字母列表的文件开始。它会生成一个sed脚本(sed当然是使用 ),在行尾连续添加一个氨基酸字母,打印它,删除它,然后移动到下一个字母。

肽-A.sh

printf "%s\n" 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.5Bpeptides-B.sh在我的磁盘上调用了脚本),并检查了它peptides.5A并且peptides.5B是相同的。

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


编辑行首而不是行尾会产生大约 20% 的性能改进。

代码:

printf "%s\n" 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 位数字时,原版开始变得不稳定。

于 2015-10-16T22:47:29.430 回答
1

免责声明:虽然我很高兴能根据 base-19 数字计算出这个算法,但它的速度非常慢(3 个字母的字符串需要 8 秒,4 个字母的字符串需要 160 秒,都包含 19 个氨基酸,运行在与 Jonathan Leffler 暗示的其他解决方案相比,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)。

然后我们从Nto循环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
于 2015-10-16T21:58:04.423 回答
1

crunch可用于 Kali 发行版

crunch 7 7 ADEFGHIKLMNPQRSTVWY
于 2018-06-08T15:37:27.573 回答