【生物信息学】嵌套循环的一些问题，打赏大神^.^

rougayo

首先介绍下我为什么要写这个脚本：

我在LINUX上运行一个名为prosplign的比对软件（可以通过蛋白序列与基因组序列比较来找出基因组序列上的exon区域），它的运行代码为：prosplign -full -nfa genome_file(一条基因组序列) -pfa  protein_file（一条蛋白序列） -o result_file(比对结果输出的文件夹)   这个程序一次只能运行一条基因组序列及一条蛋白序列，所以比较麻烦，才想着写个脚本来解决。
genome_file 在一个文件夹内，有300多条（如图一）。protein_file在另一个文件夹内，有60多条（如图二）。

我的目的是能够将genome_file中的文件与protein_file中的文件两两匹配，最后能自动用于prosplign的运行代码中，从而获得目标结果。

rougayo

喔对 我的思路是这样的。
genome和protein的序列信息之前都是分别在一个文本文件中，于是我分别将它们分割了，每个序列独立一个文件，组成一个文件夹。
然后写了一个嵌套循环的脚本，想要将它们两两结合。
最后再写一个脚本来运行结合后的genome、protein，来运行prosplign。
附上我写的嵌套循环脚本（能运行，但是输出的文件是空的，T.T）

1 #!/usr/bin/perl
      2 use strict;
      3 use warnings;
      4
      5 my$nfainput=shift||die;
      6 my$pfainput=shift||die;
      7 my$output=shift||die;
      8
      9 open NFA,"<",$nfainput||die"$!\n";
     10 open PFA,"<",$pfainput||die"$!\n";
     11 my $i=1;
     12 my @nfaa=<NFA>;
     13 my @pfaa=<PFA>;
     14 foreach my $nfa(@nfaa){
     15   foreach my $pfa(@pfaa){
     16 open  OUTPUT,">",$output.'/'.$i;
     17 print OUTPUT$_;
     18 close OUTPUT;
     19 $i++;
     20 }
     21 }
     22 close (NFA);
     23 close (PFA);

jason680

回复 1# rougayo

$ ls gnl_*.fa
gnl_1.fa  gnl_2.fa  gnl_3.fa  gnl_4.fa  gnl_5.fa

$ ls ENS*_*.fa
ENSALP00000050043_4.fa  ENSFCAP0000009652_1.fa
ENSAMEP0000007250_3.fa  ENSGGOP0000007456_2.fa



$ ls gnl_*.fa | awk -vrun=0 -F_ 'BEGIN{cmd="ls ENS*_*.fa";while(cmd | getline)a[$NF]=$0;close(cmd)}{if(!a[$NF]){print "cannot find ENS*_"$NF" file\n";next}cmd="prosplign -full -nfa "$0" -pfa "a[$NF]" -o result_"$NF".txt";print "cmd="cmd;if(run)system(cmd)}'
cmd=prosplign -full -nfa gnl_1.fa -pfa ENSFCAP0000009652_1.fa -o result_1.fa.txt
cmd=prosplign -full -nfa gnl_2.fa -pfa ENSGGOP0000007456_2.fa -o result_2.fa.txt
cmd=prosplign -full -nfa gnl_3.fa -pfa ENSAMEP0000007250_3.fa -o result_3.fa.txt
cmd=prosplign -full -nfa gnl_4.fa -pfa ENSALP00000050043_4.fa -o result_4.fa.txt
cannot find ENS*_5.fa file

rougayo

cannot find ENS*_77.fa file

cannot find ENS*_78.fa file

cannot find ENS*_79.fa file

cmd=prosplign -full -nfa gnl_7.fa -pfa ENSMODP00000020757_7.fa -o result_7.fa.tx   t
cannot find ENS*_80.fa file

cannot find ENS*_81.fa file

cannot find ENS*_82.fa file


大神您好！我将gnl_*.fa替换为bos-genome这个文件夹，ENS*_*.fa替换为HBA1_fasta这个文件夹。

报错了“zsh: command not found: $” 。于是我便将$给删除了。
最后运行出来的结果如上。
我想请教一下，因为基因组序列有300多条 而蛋白序列只有60多条。我试了您的方法 ，发现基因组序列只能比对上60多条，就是多少条基因组序列匹配多少条蛋白序列。我想问有什么方法能使每一条基因组序列都能与那66条蛋白序列想匹配一次呢？就是会出现300*66种匹配结果。

moperyblue

1. for f_genome in $(ls genome); do
   
2.     f1=${f_genome%.*}
   
3.     for f_protein in $(ls protein); do
   
4.         f2=${f_protein%.*}
   
5.         prosplign -full -nfa "$f_genome" -pfa  "$f_protein" -o result_"${f1}_$f2".txt
   
6.     done
   
7. done

复制代码