筛选具有与另一个文件中某个字符串匹配的序列的fasta文件

Question

使用BLAST，我获得了一个文件，其中有两列用制表符分隔，一列有物种名称，另一列有基因名称(参考数据库中最相似基因的名称)。我的目标是在第一个文件中找到相关基因名称最常见的所有物种名称，并使用这个物种列表来过滤第二个文件(FASTA格式)，只保留这些物种。

例如，如果这是BLAST文件：

sp1 XXX
sp2 AAA
sp3 XXX
sp4 XXX
sp5 BBB

与最常见基因匹配的物种将是sp1、sp3和sp4。然后在最初包含许多物种的FASTA文件中，我只想保留物种sp1、sp3和sp4的序列。从这个开始：

>sp1
AATCGAGTCGT
>sp2
AATCGCGTCGT
>sp3
AATCGAGTAGT
>sp4
AATCGAGTCGG
>sp5
AAACGAGTCGT

到此fasta文件：

>sp1
AATCGAGTCGT
>sp3
AATCGAGTAGT
>sp4
AATCGAGTCGG

到目前为止，我已经尝试了几天不同的方法，但出于某种原因，它永远不会奏效。我尝试创建这个脚本(don.awk)来获得最常见的基因名称(第二列)：

BEGIN {
FS="t"
}
{a[$2]++; if (a[$2] > comV) { comN=$2; comV=a[$2]} }
END {
printf("%sn", comN, comV)
}

然后用运行

nawk -f don.awk blastfile

然后我尝试将名称分配给一个变量(例如var1)，并在中使用它

awk -F't' '$2 ~ /$var1/ {print $0}' blastfile > result

首先过滤原始文件。但由于某些原因，我无法保存变量，或者出现错误。但我想肯定有更简单的方法。

Answer 1

使用GNUawk和grep:

grep -A 1 -f <(awk '
BEGIN{PROCINFO["sorted_in"] = "@val_num_desc"}
{a[$2]++; b[$1]=$2}
END{for (i in a){j++; if (j==1){tgt=i; break}}
for (ii in b){if(tgt==b[ii]){print ">"ii}}}' blast.txt) input.fasta | grep -v '^--' > filtered.fasta

输出

>sp1
AATCGAGTCGT
>sp3
AATCGAGTAGT
>sp4
AATCGAGTCGG

描述

• 将数组的排序顺序设置为降序，数组值的数值：

BEGIN{PROCINFO["sorted_in"] = "@val_num_desc"}

• 使用值作为索引填充第2列值的数组计数。并使用第1列作为关键字，第2列作为值填充第二个数组

{a[$2]++; b[$1]=$2}

• 提取第一个数组索引作为目标基因名称，用于从第二个数组中筛选物种名称

for (i in a){j++; if (j==1){tgt=i; break}}

• 使用上面提取的基因名称获取物种名称，用于过滤输入。asta文件：

for (ii in b){if(tgt==b[ii]){print ">"ii}}

• 使用grep、-A 1开关和-f开关提取所需的fasta记录

grep -A 1 -f <(awk...) input.fasta`

• 将|管道传输到另一个grep，以消除连续匹配组之间的--，并将输出重定向到文件

| grep -v '^--' > filtered.fasta

---

-或-仅使用GNUawk:

awk 'BEGIN{PROCINFO["sorted_in"] = "@val_num_desc"}
BEGINFILE{fnum++}
fnum==1{a[$2]++; b[$1]=$2}
ENDFILE{if(fnum==1){
for (i in a){j++; if (j==1){tgt=i; break}}} 
for (ii in b){if(tgt!=b[ii]){delete b[ii]}}
}
fnum==2{for (i in b){if($0==">"i){
getline n; 
printf "%sn%sn", $0, n}}
}' blast.txt input.fasta

详细信息

• 将数组的排序顺序设置为降序，数组值的数值：

BEGIN{PROCINFO["sorted_in"] = "@val_num_desc"}

• 对于处理的每个文件，递增变量"fnum">

BEGINFILE{fnum++}

• 对于处理的第一个文件，使用值作为索引填充第2列值的数组计数。并使用第1列作为关键字，第2列作为值填充第二个数组

fnum==1{a[$2]++; b[$1]=$2}

• 处理完第一个文件后，提取第一个数组索引作为目标基因名称，用于从第二个数组中筛选物种名称。然后通过删除与目标基因名称不匹配的条目来过滤第二个数组

ENDFILE{if(fnum==1){
for (i in a){j++; if (j==1){tgt=i; break}}} 
for (ii in b){if(tgt!=b[ii]){delete b[ii]}}
}

• 如果处理第二个文件，则将当前行与筛选数组中的每个条目进行比较，如果匹配，则打印当前行和下一行，并将输出重定向到文件

fnum==2{for (i in b){if($0==">"i){
getline n; 
printf "%sn%sn", $0, n}}
}' blast.txt input.fasta > filtered.fasta

Answer 2

使用标准的UNIX工具加上perl一行程序来提取最频繁的基因。使用grep和cut从爆破文件中提取物种。使用grep提取FASTA标头，然后提取具有所需物种的标头，然后仅提取序列ID。最后，要按ID从FASTA文件中提取读取，请使用seqtk subseq(将其与conda一起安装)。

cut -f2 blast.txt | sort | uniq -c | sort -k1,1n | tail -n1 | perl -lane 'print $F[1];' > most_frequent_gene.txt
grep -f most_frequent_gene.txt blast.txt | cut -f1 > species_w_most_frequent_gene.txt
grep '^>' input.fasta | grep -f species_w_most_frequent_gene.txt | grep -Po '^>KS+' > seqids_w_species_w_most_frequent_gene.txt
seqtk subseq input.fasta seqids_w_species_w_most_frequent_gene.txt > seqs_w_species_w_most_frequent_gene.fasta

本例中使用的输入和输出：

head blast.txt input.fasta most_frequent_gene.txt seqids_w_species_w_most_frequent_gene.txt seqs_w_species_w_most_frequent_gene.fasta species_w_most_frequent_gene.txt
==> blast.txt <==
sp1     XXX
sp2     AAA
sp3     XXX
sp4     XXX
sp5     BBB
==> input.fasta <==
>seq1 sp1
ACTG
>seq2 sp2
ACTGA
>seq3 sp3
ACTGAC
>seq4 sp3
ACTGAC
==> most_frequent_gene.txt <==
XXX
==> seqids_w_species_w_most_frequent_gene.txt <==
seq1
seq3
seq4
==> seqs_w_species_w_most_frequent_gene.fasta <==
>seq1 sp1
ACTG
>seq3 sp3
ACTGAC
>seq4 sp3
ACTGAC
==> species_w_most_frequent_gene.txt <==
sp1
sp3
sp4

另请参阅：

• cut
• sort
• tail
• uniq

这里，GNUgrep使用以下选项：
-P:使用Perl正则表达式
-o：只打印匹配项(每行1个匹配项)，不打印整行
-f file：从文件中获取模式，每行一个。

CCD_ 22：使正则表达式引擎"；保持"；它在CCD_ 23之前匹配的所有内容，并且不将其包括在匹配中。特别是，在打印匹配项时忽略正则表达式的前一部分。

Perl单行使用以下命令行标志：
-e:neneneba告诉Perl在线查找代码，而不是在文件中查找代码
-n：一次循环输入一行，默认情况下将其分配给$_
-l：在线执行代码之前，剥去输入行分隔符(默认为*NIX上的"n")，并在打印时附加
-a:在空白处或在-F选项中指定的regex上将$_拆分为数组@F。

• perldoc perlrun：如何执行Perl解释器：命令行开关
• perldoc perlrequick:Perl正则表达式快速入门

要安装seqtk和此处未列出的许多其他生物信息学工具，请使用conda，特别是miniconda，例如：

conda create --channel bioconda --name seqtk seqtk
conda activate seqtk
# ... use seqtk here ...
conda deactivate

• conda
• conda create

Answer 3

使用任何awk:

$ cat tst.awk
NR == FNR {
gene2species[$2] = ($2 in gene2species ? gene2species[$2] FS : "") $1
if ( ++cnt[$2] > maxCnt ) {
maxCnt = cnt[$2]
maxGene = $2
}
next
}
FNR == 1 {
split(gene2species[maxGene],tmp)
for ( i in tmp ) {
maxSpecies[">"tmp[i]]
}
}
/^>/ {
isMaxSpecies = ($1 in maxSpecies ? 1 : 0)
}
isMaxSpecies

$ awk -f tst.awk blastfile fastafile
>sp1
AATCGAGTCGT
>sp3
AATCGAGTAGT
>sp4
AATCGAGTCGG

以上假设只有1个物种具有最大数量，如您提供的示例中所示。如果可以有多个物种具有相同的最大基因数，那么将其调整为：

$ cat tst.awk
NR == FNR {
gene2species[$2] = ($2 in gene2species ? gene2species[$2] FS : "") $1
if ( ++cnt[$2] > maxCnt ) {
maxCnt = cnt[$2]
}
next
}
FNR == 1 {
for ( gene in cnt ) {
if ( cnt[gene] == maxCnt ) {
split(gene2species[gene],tmp)
for ( i in tmp ) {
maxSpecies[">"tmp[i]]
}
}
}
}
/^>/ {
isMaxSpecies = ($1 in maxSpecies ? 1 : 0)
}
isMaxSpecies

例如：

$ cat blastfile
sp1 XXX
sp2 AAA
sp3 XXX
sp4 AAA
sp5 BBB

$ awk -f tst.awk blastfile fastafile
>sp1
AATCGAGTCGT
>sp2
AATCGCGTCGT
>sp3
AATCGAGTAGT
>sp4
AATCGAGTCGG

Answer 4

如果您想使用awk以外的其他工具，可以执行以下操作：

#!/bin/bash
mostCommont=$(awk '{ print $2 }' ./input | sort | uniq -c | sort -nr | head -1 | awk '{print $2}')
grep "$mostCommont" ./input | awk '{print $1}'

或者作为一行：

grep "$(awk '{ print $2 }' ./input | sort | uniq -c | sort -nr | head -1 | awk '{print $2}')" ./input | awk '{print $1}'

只需将./input替换为输入文件的路径即可。

根据您的样本数据，这将打印以下内容：

sp1
sp3
sp4

如果你想保留基因名称，只需删除最终的... | awk '{print $1}