APDB文件包含蛋白质构象的许多段落。
每个构象都以关键字开头原子并以关键字结尾结尾。
我试图在 bash 中读取文件,以便读取从 ATOM 到 END 的每一行,但我不想读取 END 这个词。
我想对每个构象(段落)执行此操作并将每个段落存储在数组中。
该文件看起来像这样:
ATOM line 1...
ATOM line 2...
ATOM line 3...
# More lines....
END
ATOM line 1...
ATOM line 2...
ATOM line 3...
# more lines...
END
一个原子到结尾是一种构象。
我希望能够将每个构象读入一个数组,包括 ATOM 但不包括 END。
我可以读取两个关键字之间的文本(不包括这两个单词),但我不知道如何包含起始词,但排除结束词。
还将每个构象读入数组,这样conf[0]
= 第一个构象 , conf[1]
= 第二个构象 依此类推不起作用。
代码:
#!/bin/bash
filename='coor.pdb'
echo Start
i=0
while read line; do
conf[$i]=$(sed -n '/ATOM/,/END/{//!p}')
i=i+1
done < $filename
echo $conf[0] > first_frame.data
答案1
#!/bin/bash
filename='coor.pdb'
echo Start
i=1
input=false
while read -r line
do
if [ "${line%% *}" == "ATOM" ]
then
input=true
elif [ "${line%% *}" == "END" ]
then
((i++)) # increase variable i by 1 == (i+1)
rm -f "${i}_frame.data" # remove ${i}_frame.data" if exist
input=false # stop output lines until next ATOM
fi
if $input # if var INPUT is true add line to ${i}_frame.data file
then
echo "$line" >> "${i}_frame.data"
fi
done < "$filename"
对于未来的一些sed提示:
sed '/ATOM/,/END/!d;/END/d'
sed -n '/ATOM/{:;N;s/\nEND//;T;p}'
所以你可以做任务:
nl -s'.frame.data' -b p"^END" coor.pdb |
sed -n '/ATOM/{s/^/echo \"/;:;s/ \{6,\}//;N;s/END//;T;s/\n */\">/p}' |
bash
答案2
bash 中的文本处理速度很慢。纯 bash 字符串操作对于变量中已有的文本或读取非常小的文件很有用。我怀疑计算生物学文件通常不会很小,因此使用这样的工具awk
启动成本很小,但处理文本的速度比 bash 快得多。
假设您真的只想拆分文件pdb
:
awk -v RS='\nEND\n' '{ fn="frame" NR ".pdb"; print > fn; close(fn) }' "$filename"
让 awk 用作\nEND\n
输入记录分隔符,而不是换行符,那么您甚至可以使用它的记录计数器。输出记录分隔符仍然是默认的ORS="\n"
。 (Costas 提出了非常好的建议。我对其进行了调整,因此END
必须位于行的开头,并添加close
以确保它不会在具有很多构象的输入上使用大量文件描述符。)
我最初的想法是:
awk 'BEGIN{i=0; fn="frame0.pdb"}
!/^END/ { print > fn; }
/^END/{ close(fn); fn="frame" ++i ".pdb"; }' \
"$filename"
awk 缓存文件句柄,因此多个print > fn
不会导致文件关闭重新打开。 (close(fn)
就是这样做的。它的存在只是为了提高效率,因此 awk 最终不会打开大量文件。)
逻辑是:将每一整行打印到当前文件名。当您看到一行时END
,请转到下一个文件名。如果在最后一行之后没有另一行END
,则永远不会写入新文件名,并且不会创建残留的最后一个文件。
OTOH,如果您想对内存中的行块数组执行某些操作:
# add a `!/^END/` condition to the concat block if you want to avoid a stray newline after each END
awk 'BEGIN{i=0}
!/^END/ { arr[i] = arr[i] $0 "\n"; } # concat onto this array element
/^END/ { i++; }
END{for (k in arr) { print arr[k]; > ("frame" k ".pdb") } }' \
"$filename"
然后你就可以在块中随意处理 awk 行数组了END
。它具有良好的正则表达式功能。
bash 驱动 sed 的尝试失败(nvm,失败是因为sed
不像 shell 那样一次读取一个字节read
):
i=0
while true; do
outf="frame${i}.data";
##### DON'T USE THIS, sed READS TOO MUCH #####
strace -o sed.tr sed '/^END/q42' > "$outf"; # strace to see that the 2nd sed invocation finds the file empty
ret=$?;
((i++));
if [[ $ret == 0 ]];then # sed didn't see END before EOF
[[ -s $outf ]] || rm -f "$outf"; # clean up empty last file
break;
elif [[ $ret != 42 ]]; then
echo some other sed error;
break;
fi;
done < "$filename"