在aaa大于300的fasta文件中选择序列,“ C”出现至少4次
收藏

我有一个包含蛋白质序列的fasta文件。我想选择超过300个氨基酸和半胱氨酸(C)氨基酸出现4次以上的序列。
我用这个命令选择了超过300aa的序列:

 cat 72hDOWN-fasta.fasta | bioawk -c fastx 'length($seq) > 300{ print ">"$name; print $seq }' 

一些序列示例:
  >jgi|Triasp1|216614|CE216613_3477
 MPSLYLTSALGLLSLLPAAQAGWNPNSKDNIVVYWGQDAGSIGQNRLSYYCENAPDVDVI
 NISFLVGITDLNLNLANVGNNCTAFAQDPNLLDCPQVAADIVECQQTYGKTIMMSLFGST
 YTESGFSSSSTAVSAAQEIWAMFGPVQSGNSTPRPFGNAVIDGFDFDLEDPIENNMEPFA
 AELRSLTSAATSKKFYLSAAPQCVYPDASDESFLQGEVAFDWLNIQFYNNGCGTSYYPSG
 YNYATWDNWAKTVSANPNTKLLVGTPASVHAVNFANYFPTNDQLAGAISSSKSYDSFAGV
 MLWDMAQLFGNPGYLDLIVADLGGASTPPPPASTTLSTVTRSSTASTGPTSPPPSGGSVP
 QWGQCGGQGYTGPTQCQSPYTCVVESQWWSSCQ* 


最佳答案:

我不知道bioawk,但我假设它与相同,有一些初始解析和常量定义。
我将按以下步骤进行。假设您希望查找长度大于300且字母Cin的4倍以上的字符串,则可以执行以下操作:

bioawk -c fastx '
   (length($seq) > 300) && (gsub("C","C",$seq)>4) {
       print ">"$name; print $seq
   }' 72hDOWN-fasta.fasta

但这假设seq是完整的字符序列。
其背后的想法如下。gsub命令在字符串中执行替换并返回它所做的全部替换。因此,如果我们用“c”替换所有字符“c”,实际上我们并没有改变字符串,而是得到字符串中“c”的总数。
POSIX standard IEEE Std 1003.1-2017开始:
gsub(ere, repl[, in]):行为类似于sub(见下文),但它应替换正则表达式的所有出现(类似于
ed实用程序全局替换)in$0或in参数中,
如果指定。
sub(ere, repl[, in ]):替换stringrepl中扩展正则表达式的第一个实例
并返回替换的数目。与号(ere
)字符串中出现的字符串in应替换为&中的字符串
这和ERE吻合。前面有一个
应解释为
字符。连续发生两次
字符应解释为单个字符
文字字符。任何其他事件
<反斜杠>(例如,在任何其他字符之前)应
被视为一个文本字符。注意,如果repl
是字符串文字(词汇标记字符串;参阅Grammar
在任何词法之后处理'abpult>字符。
处理,包括任何词法<反斜杠>转义序列
处理。如果指定了in且不是左值(请参见
Expressions in awk),行为未定义。如果省略repl,awk
应使用当前记录(in)代替。
注:BioAwk基于Brian Kernighan's awk中的"The AWK Programming Language", by Al Aho, Brian Kernighan, and Peter Weinberger (Addison-Wesley, 1988, ISBN 0-201-07981-X) 。我不确定这个版本是否与POSIX兼容。

公众号