在fasta文件中选择序列超过300 aa,"C"至少出现4次

Question

我有一个包含蛋白质序列的fasta文件.我想选择超过300个氨基酸的序列,半胱氨酸(C)氨基酸出现超过4次.

我用这个命令来选择超过300 aa的序列:

 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* 

Answer 1

我不知道,bioawk但我认为它与awk相同,具有一些初始解析和常量定义.

我将按如下方式进行.假设你想要找到字母C大于4倍且长度超过300 的字符串,那么你可以这样做:

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

但这假定这seq是完整的字符序列.

其背后的想法如下.该gsub命令在字符串中执行替换并返回它所做的总替换.因此,如果我们用"C"替换所有字符"C",我们实际上并没有改变字符串,而是返回字符串中"C"的总量.

从POSIX标准IEEE Std 1003.1-2017:

gsub(ere, repl[, in]):表现得像sub(见下文),除了它应该在指定时ed替换$0in参数中或in中的所有正则表达式(如实用程序全局替换).

sub(ere, repl[, in ]):替换字符串repl代替扩展正则表达式的第一个实例的ere字符串in ,并返回替换的数量.& 出现在字符串中的<ampersand>()repl应替换为in 与ERE匹配的字符串.前面带有<反斜杠>的<&符号>应解释为文字<&符号>字符.出现两个连续的<反斜杠>字符应解释为单个文字<反斜杠>字符.任何其他出现的<反斜杠>(例如,在任何其他字符之前)都应被视为文字<反斜杠>字符.请注意,如果repl 是字符串文字(词法标记STRING;请参阅语法),则在任何词法处理之后都会处理<&符号>字符,包括任何词法<反斜杠> - 序列处理.如果in已指定且它不是左值(请参阅 awk中的表达式),则行为未定义.如果in省略,awk将使用当前记录($0)代替它.

注意: BioAwk基于Brian Kernighan的awk,其由Al Aho,Brian Kernighan和Peter Weinberger(Addison-Wesley,1988,ISBN 0-201-07981-X)的"The AWK Programming Language"中记载 .我不确定这个版本是否与POSIX兼容.