过滤 .fasta 文件

Question

我有以下两个文件：

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

>A0A2K5RDW4|unreviewed|Olfactory receptor|taxID:2715852
MAHINCTQVAEFILVGLTDHQELKMPLFVLFLSIYLFTVVGNLGLILLIRADSRLSTPMYFFLSNLAFVDFCYSSVITPK
MLGTFLYKQNVISFNACATQLCCFLTFMISECLLLASIAYDRYVAICNPLLYMVVMSPGICIQLVAVPYSYSFLMALFHT
ILTFRLSYCHSNIVNHFYCDDMPLLRLTCSDTSSKQLWIFACAGIMFTSSLLIVFVSYMFIISAILRMRSAEGRQKAFST
CGSHMLAVTIFYGTLIFMYLQPSSNHALDTDKMASVFYTVIIPMLNPLIYSLRNKEVKEALKKVIINKNQ
>A0A2K5RFA5|unreviewed|Olfactory receptor|taxID:2715852
MSRWSNGSSNVSYTSFLLAGFPGLRESRALLVLPFLSLYLMIIFTNALVIHTVLTQQSLHQPMYLLIALLLAVNICATST
VVPAMLFSFSTHFNCISLARCLIQMFCIYFLIVFDCNILLVMALDRYVAICYPLCYPEIVTGQLLAGLVGAAATRSTCIV
APVVFLASRVHFCRSDMIQHFACEHMALMKLSCGDISLNKTVGLTVRIFNRVLDMLLLGTSYSHIIHAAFRISSGRARSK
ALNTCGSHLLVIFTVYLSTMSSSIVYRTAHTASQDVHNLLSTFYLLLPCLVNPIIYGTRTKEIRYHLGEQFLSSGPRDTM
MSI

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

A0A2K5BUY4
A0A2K5BUY7
A0A2K5BUZ1
A0A2K5BW26
A0A2K5BW36
A0A2K5BWA0

我想创建第二个 .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

是一个好主意，但您的代码存在几个问题：

使用
```
seq_id = substr($0, 2)
```
，您不仅可以分配 ID，还可以分配整个标题减去前导
```
>
```
。
使用
```
accessions[$1]
```
添加到数组
```
accessions
```
的索引是 ID 不带前导
```
>
```
。但
```
seq_id
```
```
 中 
```
seq_id in accessions 的值具有领先
```
>
```
。
您决定在处理标题时进行打印（在
```
if (substr($0, 1, 1) == ">" {...}
```
内），因此序列的其余部分尚未在您的
```
seq
```
变量中。
...

如果标题行始终以

开头，并且没有其他行以

开头，您可以尝试：

awk -F '|' '
  NR == FNR {accessions[">" $1]; next}
  /^>/ {to_print = $1 in accessions}
  to_print' platyrrhini_or_accessions.txt platyrrhini_784ORs.fasta

awk

语句是应用于输入记录的模式-操作对，当且仅当当前记录的模式为真时才执行该操作。在您的情况下，记录是常规线路。

第一个语句 (
```
NR == FNR
```
) 的模式仅适用于第一个文件的行。它的操作 (
```
accessions[">" $1]; next
```
) 存储在
```
platyrrhini_or_accessions.txt
```
中找到的 ID，并添加前导
```
>
```
作为
```
accessions
```
数组的索引，并移动到下一行 (
```
next
```
)，从而跳过其他 2 个语句.
第二条语句 (
```
/^>/
```
) 的模式是一个正则表达式，对于所有以
```
>
```
开头的行，其计算结果为 true。如果该行的第一个字段（即以
```
to_print = $1 in accessions
```
开头的 ID）是
```
to_print
```
数组的索引，则操作 (
```
>
```
) 将 true 分配给变量
```
accessions
```
，否则为 false。
第三条语句只有模式部分（
```
to_print
```
），没有动作部分。因此，如果
```
to_print
```
为 true，即如果我们遇到的最后一个 ID 位于
```
platyrrhini_or_accessions.txt
```
中，则它将默认操作（打印）应用于当前行。

注意：

awk

中没有真正的布尔类型； 0 和空字符串被视为 false，其他任何内容都被视为 true。对于任何符合 POSIX 标准的

awk

，在上面的代码中，

to_print

变量采用值 0（假）或非零数值（真），可能为 1（但这不是 POSIX 指定的）。

过滤 .fasta 文件

问题描述投票：0回答：1

1个回答

最新问题

过滤 .fasta 文件

问题描述 投票：0回答：1

1个回答

最新问题

问题描述投票：0回答：1