过滤 JSON 文件并使用输出来过滤 .fasta 文件

Question

我有以下两个文件：

1. 包含 >700 个蛋白质序列的 .fasta 文件。每个序列占用几行，并有一个标题，其中包含有关蛋白质的一些信息，包括登录代码（例如 A0A2K5BUY4）。一个完整的序列（包括标题）如下所示：

   >A0A2K5BUY4|unreviewed|Olfactory receptor|taxID:37293
   MKNRSMEIEFILLGLTDNPQLQIVIFLFLFLNYTLSLMGNLVIIILTLLDPRLKTPMYFFLRNFSFLEIIFTTVCIPRFL
   ITIVTRDKTISYNNCVTQLFFILLPGVTEFYLLAAMSYDRYVAICKPLRYPIIMSSKVCYQLVLSSWVTGFLIIFPPLVM
   GLKLDFCASKTIDHFMCETSPILQISCTDTHVLELMSFILAVVTLVITLVLVILSYTCIIKTILKFSSAQQRSKAFSTCT
   SHMIVVSMTYGSCIFMYIKPSAKERVTVSKGVALLYTSIAPLLNPFIYTLRNQQVKEVFWDILQKNLCFSKRPF
   

2. 一个 .txt 文件，有 229 行，每行包含一个登录代码，就像 .fasta 文件中序列标题中的代码一样。

我想创建第二个 .fasta 文件，该文件仅包含其登录码位于 .txt 文件中的序列（包括标题）。

我对生物信息学和一般编程都很陌生，所以我在 chatGPT 的帮助下做了以下工作：

# Use pcregrep to extract sequences that match the accession codes
pcregrep -M -f platyrrhini_or_accessions.txt platyrrhini_784ORs.fasta | \
  pcregrep -M ">([^>]+)(\n[^>]+)*" > filtered_platyrrhini_ORs.fasta

但这并没有捕获序列，因为输出仅包含匹配蛋白质的标题：

>A0A2R8PLQ0|unreviewed|Olfactory receptor|taxID:9483
>A0A5F4VRG8|unreviewed|Olfactory receptor|taxID:9483
>A0A5F4VRT9|unreviewed|Olfactory receptor|taxID:9483
>A0A5F4VS80|unreviewed|Olfactory receptor|taxID:9483

我也尝试过使用 grep，但它不起作用，因为序列不占用固定的行数：

grep -A 1 -f platyrrhini_or_accessions.txt platyrrhini_784ORs.fasta > filtered_platyrrhini_ORs.fasta

这捕获了 229 个标头，但每个序列仅捕获 1 行。

最后我使用了 awk：

awk 'BEGIN {FS="|"} NR==FNR {accessions[$1]; next} 
     {if (substr($0, 1, 1) == ">") {seq_id = substr($0, 2); seq = $0} 
      else {seq = seq "\n" $0} 
      if (seq_id in accessions) {print seq_id "\n" seq}}' \
platyrrhini_or_accessions.txt platyrrhini_784ORs.fasta > filtered_platyrrhini_ORs.fasta

输出文件为空。

我真的很感激任何帮助。

Answer 1

如果标题行始终以

>

开头，并且没有其他行以

>

开头，您可以尝试：

awk -F '|' 'NR == FNR {accessions[">" $1]; next}
  /^[>]/ {ok = ($1 in accessions) ? 1 : 0} ok' platyrrhini_or_accessions.txt platyrrhini_784ORs.fasta