我希望生成一个文件夹,包含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.5A
和peptides.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非常正确地评论了输出大小)。
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