使用 awk 打印标头名称和子字符串

Question

>我尝试使用此代码打印基因名称的标头，然后根据其位置提取子字符串，但它不起作用

>output_file
cat input_file | while read row; do
echo $row > temp
geneName=`awk '{print $1}' tmp`
startPos=`awk '{print $2}' tmp`
endPOs=`awk '{print $3}' tmp`
for i in temp; do
echo ">${geneName}" >> genes_fasta ;
echo "awk '{val=substr($0,${startPos},${endPOs});print val}' fasta" >> genes_fasta
done
done

input_file

nad5_exon1 250405 250551
nad5_exon2 251490 251884
nad5_exon3 195620 195641
nad5_exon4 154254 155469
nad5_exon5 156319 156548

法斯塔

atgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgc............

这是我错误的输出文件

>
awk '{val=substr(pull_genes.sh,,);print val}' unwraped_carm_mt.fasta
>
awk '{val=substr(pull_genes.sh,,);print val}' unwraped_carm_mt.fasta
>
awk '{val=substr(pull_genes.sh,,);print val}' unwraped_carm_mt.fasta
>
awk '{val=substr(pull_genes.sh,,);print val}' unwraped_carm_mt.fasta
>
awk '{val=substr(pull_genes.sh,,);print val}' unwraped_carm_mt.fasta
>
awk '{val=substr(pull_genes.sh,,);print val}' unwraped_carm_mt.fasta

输出应如下所示：

>name1
atgcatgcatgcatgcatgcat
>name2
tgcatgcatgcatgcat
>name3
gcatgcatgcatgcatgcat
>namen....

Answer 1

您可以通过对awk的单次调用来做到这一点，这将比在 shell 脚本中循环并每次迭代调用 4 次awk效率高几个数量级。由于您有 bash，您可以简单地使用命令替换并将fasta的内容重定向到awk变量，然后简单地从fasta文件中输出包含开头到结束字符的标题和子字符串。

例如：

awk -v fasta=$(<fasta) '{print ">" $1; print substr(fasta,$2,$3-$2+1)}' input

或在BEGIN规则中使用getline：

awk 'BEGIN{getline fasta<"fasta"}
{print ">" $1; print substr(fasta,$2,$3-$2+1)}' input

示例输入文件

注意：开始值和结束值已减少，以适应示例的 129 个字符：

$ cat input
rad5_exon1 1 17
rad5_exon2 23 51
rad5_exon3 110 127
rad5_exon4 38 62
rad5_exon5 59 79

示例的前 129 个字符fasta

$ cat fasta
atgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgcatgc

示例使用/输出

$ awk -v fasta=$(<fasta) '{print ">" $1; print substr(fasta,$2,$3-$2+1)}' input
>rad5_exon1
atgcatgcatgcatgca
>rad5_exon2
gcatgcatgcatgcatgcatgcatgcatg
>rad5_exon3
tgcatgcatgcatgcatg
>rad5_exon4
tgcatgcatgcatgcatgcatgcat
>rad5_exon5
gcatgcatgcatgcatgcatg

仔细查看，让我知道我是否理解您的问题要求。如果您对解决方案还有其他疑问，也请告诉我。

Answer 2

如果我理解正确，怎么样：

awk 'NR==FNR {fasta = fasta $0; next}
{
printf(">%s %sn", $1, substr(fasta, $2, $3 - $2 + 1))
}' fasta input_file > genes_fasta

• 它首先读取fasta文件并将序列存储在变量fasta中。
• 然后它逐行读取input_file，提取从$2开始且长度为$3 - $2 + 1的fasta子字符串。(请注意，substr函数的第 3 个参数是长度，而不是 endpos。

希望这有帮助。

Answer 3

让它工作了！ 这是用于从 fasta 文件中提取子字符串的脚本

cat genes_and_bounderies1 | while read row; do
echo $row > temp
geneName=`awk '{print $1}' temp`
startPos=`awk '{print $2}' temp`
endPos=`awk '{print $3}' temp`
length=$(expr $endPos - $startPos)
for i in temp; do
echo ">${geneName}" >> genes_fasta
awk -v S=$startPos -v L=$length '{print substr($0,S,L)}' unwraped_${fasta} >> genes_fasta
done
done