忘记密码   免费注册 查看新帖 | 论坛精华区

ChinaUnix.net

  平台 论坛 博客 认证专区 大话IT 视频 徽章 文库 沙龙 自测 下载 频道自动化运维 虚拟化 储存备份 C/C++ PHP MySQL 嵌入式 Linux系统
最近访问板块 发新帖
查看: 502 | 回复: 4

【生物信息学】嵌套循环的一些问题,打赏大神^.^ [复制链接]

论坛徽章:
0
发表于 2017-02-28 14:03 |显示全部楼层
首先介绍下我为什么要写这个脚本:

我在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的运行代码中,从而获得目标结果。



图一

图一

图二

图二
打赏鼓励一下!

论坛徽章:
0
发表于 2017-02-28 14:19 |显示全部楼层
喔对 我的思路是这样的。
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);

论坛徽章:
109
技术图书徽章
日期:2013-10-01 15:32:13戌狗
日期:2013-10-25 13:31:35金牛座
日期:2013-11-04 16:22:07巳蛇
日期:2014-05-09 16:43:18巨蟹座
日期:2014-10-23 17:48:38子鼠
日期:2013-11-18 18:48:57白羊座
日期:2013-11-29 10:09:11狮子座
日期:2013-12-12 09:57:42白羊座
日期:2013-12-24 16:24:46辰龙
日期:2014-01-08 15:26:12技术图书徽章
日期:2014-01-17 13:24:40巳蛇
日期:2014-02-18 14:32:59
发表于 2017-02-28 14:29 |显示全部楼层
回复 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

论坛徽章:
0
发表于 2017-02-28 15:00 |显示全部楼层
本帖最后由 rougayo 于 2017-02-28 15:39 编辑


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种匹配结果。




论坛徽章:
26
15-16赛季CBA联赛之八一
日期:2016-02-22 19:10:4215-16赛季CBA联赛之青岛
日期:2016-11-26 17:00:4615-16赛季CBA联赛之深圳
日期:2016-12-01 10:34:0415-16赛季CBA联赛之新疆
日期:2016-12-07 10:24:2915-16赛季CBA联赛之同曦
日期:2016-12-15 12:06:43CU十四周年纪念徽章
日期:2016-12-18 13:03:4415-16赛季CBA联赛之吉林
日期:2017-01-03 15:52:2515-16赛季CBA联赛之辽宁
日期:2017-01-04 14:58:2415-16赛季CBA联赛之辽宁
日期:2017-01-15 09:42:512016科比退役纪念章
日期:2017-02-06 17:21:50黑曼巴
日期:2017-02-10 15:46:1215-16赛季CBA联赛之上海
日期:2017-03-18 10:14:54
发表于 2017-02-28 19:17 |显示全部楼层

  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
复制代码
您需要登录后才可以回帖 登录 | 注册

本版积分规则

  

北京皓辰网域网络信息技术有限公司. 版权所有 京ICP证:060528号 北京市公安局海淀分局网监中心备案编号:1101082001
广播电视节目制作经营许可证(京) 字第1234号 中国互联网协会会员  联系我们:
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP