如何使用linux命令提取测序数据

Question

我想提取某些行及其后续测序数据。

有一个ecoli.ffn文件如下：

$head ecoli.ffn
>ecoli16:g027092:GCF_000460315:gi|545267691|ref|NZ_KE701669.1|:551259-572036
ATGAGCCTGATTATTGATGTTATTTCGCGT
AAAACATCCGTCAAACAAACGCTGATTAAT
>ecoli16:g000011:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGAT
>ecoli16:g000012:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGAT
CTGACAGCTGTTCTTACACTGGATTCAACC
CTGACAGCTGTTCTTACACTGGATTCAACC

和一个index.txt如下

$head index.txt
g000011
g000012

我想做的是"extract index.txt from ecoli.ffn"，理想的输出是：

>ecoli16:g000011:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGAT
>ecoli16:g000012:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGAT
CTGACAGCTGTTCTTACACTGGATTCAACC
CTGACAGCTGTTCTTACACTGGATTCAACC

我该怎么做？

Answer 1

awk 救援！

$ awk -F: -v RS=">" 'NR==FNR{n=split([=10=],t,"\n");
                             for(i=1;i<n;i++) a[t[i]]; 
                             next} 
                      in a{printf "%s", RS [=10=]}' index file  

>ecoli16:g000011:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGATCTGACAGCTGTTCTTACACTGGATTCAACC
>ecoli16:g000012:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGATCTGACAGCTGTTCTTACACTGGATTCAACC

更新请注意，这并不取决于每条记录有多少行。对于更新的输入文件，相同的脚本将为您提供此输出

$ awk -F: -v RS=">" 'NR==FNR{n=split([=11=],t,"\n");
                             for(i=1;i<n;i++) a[t[i]];
                             next}
                      in a{printf "%s", RS [=11=]}' index file

>ecoli16:g000011:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGAT
>ecoli16:g000012:55989:gi|218693476|ref|NC_011748.1|:1128430-1131042
GTGTACGCTATGGCGGGTAATTTTGCCGAT
CTGACAGCTGTTCTTACACTGGATTCAACC
CTGACAGCTGTTCTTACACTGGATTCAACC

Answer 2

使用 awk 写一个简单的脚本 ecoli.sh:

#!/bin/bash
a=`cat index.txt`
for i in $a
do
    cat ecoli.ffn|awk -F: -v i="$i" 'BEGIN{flag=0} {if( == i){print [=10=];flag=1;} if(flag ==1 &&  != i){print [=10=]; flag=0;} }'
done

那么您需要在 shell 中运行这个脚本。

Answer 3

此脚本可用于根据列表或文件的 ID 过滤 FASTA 文件，这似乎是您在这里要求的：

https://github.com/jorvis/biocode/blob/master/fasta/filter_fasta_by_ids.pl

如何使用linux命令提取测序数据

How to use linux command to extract sequencing data

linux

bioinformatics