在bash中按指定的序列长度过滤掉FASTA文件答案

【问题标题】：Filter out FASTA files by specified sequence length in bash在bash中按指定的序列长度过滤掉FASTA文件
【发布时间】：2020-05-28 13:59:55
【问题描述】：

有一个 FASTA 文件 assembly.fasta 包含 contig 名称和相应的序列：

>contig_1
CCAATACGGGCGCGCAGGCTTTCTATCGCGCGGCCGGCTTCGTCGAGGACGGGCGGCGCA
AGGATTACTACCGCAGCGGC
>contig_2
ATATAAACCTTATTCATCGTTTTCAGCCTAATTTTCCATTTAACAGGGATGATTTTCGTC
AAAATGCTGAGGCTTTACCAAGATTTTCTACCTTGCACCTTCAGAAAAAAATCATGGCAT
TTATAGACGAAATTCTCGAGAAA
>contig_3
CGTGATCTCGCCATTCGTGCCG

我只想获得长度超过 30 个字母的重叠群，并获得一个新的 FASTA 文件 assembly.filtered.fasta，其中仅包含具有重叠群名称的长序列，格式如下：

>contig_1
CCAATACGGGCGCGCAGGCTTTCTATCGCGCGGCCGGCTTCGTCGAGGACGGGCGGCGCA
AGGATTACTACCGCAGCGGC
>contig_2
ATATAAACCTTATTCATCGTTTTCAGCCTAATTTTCCATTTAACAGGGATGATTTTCGTC
AAAATGCTGAGGCTTTACCAAGATTTTCTACCTTGCACCTTCAGAAAAAAATCATGGCAT
TTATAGACGAAATTCTCGAGAAA

【问题讨论】：

能否请您添加更多详细信息如何在您的问题中实现示例输出，因为它不清楚（不是我的反对票）。
@RavinderSingh13 是的，更新了！

标签： awk fasta

【解决方案1】：

使用gnu-awk，您可以使用这个更简单的版本：

awk -v RS='>[^\n]+\n' 'length() >= 30 {printf "%s", prt $0} {prt = RT}' file

>contig_1
CCAATACGGGCGCGCAGGCTTTCTATCGCGCGGCCGGCTTCGTCGAGGACGGGCGGCGCA
AGGATTACTACCGCAGCGGC
>contig_2
ATATAAACCTTATTCATCGTTTTCAGCCTAATTTTCCATTTAACAGGGATGATTTTCGTC
AAAATGCTGAGGCTTTACCAAGATTTTCTACCTTGCACCTTCAGAAAAAAATCATGGCAT
TTATAGACGAAATTCTCGAGAAA

【讨论】：

这个解决方案输出了一些在其他两个解决方案的输出中不存在的重叠群，尽管它们是正确的，不知道
这个解决方案在某种程度上是错误的，因为如果我在此之上应用 RavinderSingh13 的解决方案，我仍然会过滤掉一些重叠群
不知道怎么会出错，这个输出与您的预期输出不匹配吗？
不，它没有。我应用了 3 种不同的解决方案，包括这个，2 个产生相同的结果，另一个产生相同的结果。我不确定哪里错了，因为真实数据太大而无法手动检查错误
这与您提供的样本数据输出相匹配。如果你能给我提供更长的样本，那么我肯定可以更深入地挖掘。

【解决方案2】：

实现您所追求的一个非常快速的方法是：

awk -v n=30 '/^>/{ if(l>n) print b; b=$0;l=0;next }
            {l+=length;b=b ORS $0}END{if(l>n) print b }' file

您可能还对BioAwk 感兴趣，它是 awk 的改编版本，用于处理 FASTA 文件

bioawk -c fastx -v '(length($seq)>30){print ">" $name ORS $seq}' file.fasta

注意： 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。

【讨论】：

如果您投了反对票，我想要求放置 cmets，否则如果有人不知道拒绝投票的原因，一个人如何改善？解决方案效果很好（对于给定的例子）并且很好地解释了然后不确定什么可能是沉默投票的原因:(

【解决方案3】：

您能否尝试使用所示示例进行跟踪、测试和编写。

awk '
/^>/{
  if(sign_val && strLen>=30){
    print sign_val ORS line
  }
  strLen=line=""
  sign_val=$0
  next
}
{
  strLen+=length($0)
  line=(line?line ORS:"")$0
}
END{
  if(sign_val && strLen>=30){
    print sign_val ORS line
  }
}
'  Input_file

说明：为上面添加详细说明。

awk '                                ##Starting awk program from here.
/^>/{                                ##Checking condition if line starts from > then do following.
  if(sign_val && strLen>=30){        ##Checking if sign_val is SET and steLen is SET then do following.
    print sign_val ORS line          ##Printing sign_val ORS and line here.
  }
  strLen=line=""                     ##Nullify variables steLen and line here.
  sign_val=$0                        ##Setting sign_val to current line here.
  next                               ##next will skip all further statements from here.
}
{
  strLen+=length($0)                 ##Checking length of line and keep adding it here.
  line=(line?line ORS:"")$0          ##Creating line variable and keep appending it to it with new line.
}
END{                                 ##Starting END block of this program from here.
  if(sign_val && strLen>=30){        ##Checking if sign_val is SET and steLen is SET then do following.
    print sign_val ORS line          ##Printing sign_val ORS and line here.
  }
}
' Input_file                         ##mentioning Input_file name here.

【讨论】：

我检查了我的数据，它的工作方式符合我的预期。