Chinaunix

标题: 求助 perl [打印本页]

作者: zhagnqiang829 时间: 2013-03-13 10:05
标题: 求助 perl
下面格式内容，把中间的大写字母的配对率测出来，
根据中间配对率大小：把小于0.6的存入文件1.txt，0.6-0.9之间的存入文件2.txt，大于0.9的存入3.txt文件

I001 195 TTGGT 205 225
E012 148 TTCGT 158 200

I001 1124 TTTCTT 1151 2978
E012 4 TAATGT 31 128

I002 17 AATTGACT 29 54
E023 5 AATTGACT 17 128

I003 16 TCGACTTAAT 35 65
E034 44 TGGACATAAT 63 200

配对率定义：设最佳匹配长度为m bp，若有k个碱基完全配对，则配对率为k/m。例如第一组序列配对率为：4/5=0.8，第二组为：2/6=0.3。。。
结果如下：
1.txt中有
I001 1124 TTTCTT 1151 2978
E012 4 TAATGT 31 128
2.txt中有：
I001 195 TTGGT 205 225
E012 148 TTCGT 158 200
I003 16 TCGACTTAAT 35 65
E034 44 TGGACATAAT 63 200
3.txt中有：
I002 17 AATTGACT 29 54
E023 5 AATTGACT 17 128

作者: rubyish 时间: 2013-03-13 10:54
本帖最后由 rubyish 于 2013-03-13 06:58 编辑

如下：

#!/usr/bin/perl
my @f = map { open my ($f), '>', $_; $f } qw/1.txt 2.txt 3.txt/;
while (<DATA>) {
next if /^\s*$/;
my $S2 = <DATA>;
my $s1 = (split)[2];
my ($s2) = $S2 =~ /([ATGC]+)/;
my ( $l, $m ) = length $s1;
for my $i ( 0 .. $l - 1 ) {
substr( $s1, $i, 1 ) eq substr( $s2, $i, 1 ) && $m++;
}
my $p = $m / $l;
my $i = $p < 0.6 ? 0 : $p < 0.9 ? 1 : 2;
print { $f[$i] } $_, $S2;
}
__DATA__
I001 195 TTGGT 205 225
E012 148 TTCGT 158 200
I001 1124 TTTCTT 1151 2978
E012 4 TAATGT 31 128
I002 17 AATTGACT 29 54
E023 5 AATTGACT 17 128
I003 16 TCGACTTAAT 35 65
E034 44 TGGACATAAT 63 200

复制代码

作者: zhagnqiang829 时间: 2013-03-13 11:54
大侠再:wink: 帮我看看，怎么产生随机序列？
下面格式内容，把中间的大写字母的GC含量测出来，
根据中间GC含量大小：把I001，I002，。。。。等开头的序列换成随机序列（此序列要和原序列GC含量和长度一样）,E012....等E开头的序列不变。
16191_NC_003284
I001CACTCA
E012ATGACAAGGAAGCGGCG

16192_NC_003284
I001CACAAGA
E012TCTACAACACGTCTGCTC
I002CATG
E023GCTGCAGAGA
GC含量定义：设序列长度为m bp，若序列有k个G+C（就是C和G的总是），则GC含量为k/m。
例如I001序列GC含量为：3/6=0.5，第二个I001序列为：3/7=0.43，第三序列为0.25,I002序列为0.5。
结果如下：
16191_NC_003284
I001AGCTTC
E012ATGACAAGGAAGCGGCG

16192_NC_003284
I001GATCAGT
E012TCTACAACACGTCTGCTC
I002TGGA
E023GCTGCAGAGA

回复 2# rubyish

作者: rubyish 时间: 2013-03-13 13:40
看看:

#!/usr/bin/perl
my @G = qw/G C/;
my @A = qw/A T/;
print /^I/ ? Rand() : $_ while <DATA>;
sub Rand {
my ( $h, $t ) = /^(.*\d+)(.*)$/;
my @t = map /[GC]/ ? $G[rand 2] : $A[rand 2], $t =~ /./g;
@t = map splice(@t, rand @t, 1), 0 .. $#t;
join '', $h, @t, $/;
}
__DATA__
16191_NC_003284
I001CACTCA
E012ATGACAAGGAAGCGGCG
16192_NC_003284
I001CACAAGA
E012TCTACAACACGTCTGCTC
I002CATG
E023GCTGCAGAGA

复制代码

作者: zzwevil 时间: 2013-03-13 14:00
新手来试试
#!/usr/bin/perl
$src='e:\perl\src1.txt';
open OPF ,$src or die "cant open ";
$opkey=0;
while (<OPF>){
if(/^$/) {next};
if ($opkey==0){
      $first=$_;
      chomp $_;
      /(\w\d+)\ +(\d+)\ +(\w+)\ +(\d+)\ +(\d+)/;
      $a=$3;
      $opkey++;
      }
elsif ($opkey==1){
      $second=$_;
      chomp $_;
      /(\w\d+)\ +(\d+)\ +(\w+)\ +(\d+)\ +(\d+)/;
      $b=$3;
      $leng=length $b;
      @aa=(split //,$a);
      @bb=(split //,$b);

      for ($i=0;$i<$leng;$i++){
            $sfta=shift @aa;
            $sftb=shift @bb;
            if ($sfta eq $sftb){$sum +=1;}
            }
      $percent=$sum/$leng;
      if($percent <=0.6){
open GO,">>","d:\\proj1\\1.txt";
         print GO $first;
                              print GO $second;
         }
      if(($percent<=0.9)&&($percent>=0.6)){
         open GO,">>","d:\\proj1\\2.txt";
                    print GO $first;
                              print GO $second;
                              }
      if($percent >=0.9){
         open GO,">>","d:\\proj1\\3.txt";
         print GO $first;
                              print GO $second;
                              }
      $opkey=0;
      }
$sum=0;
}

作者: zhagnqiang829 时间: 2013-03-13 15:44
谢谢你，写的很详细。。回复 5# zzwevil

作者: zhagnqiang829 时间: 2013-03-13 15:49
这么简练的代码就搞定了，就是要懂很费劲呀。
非常感谢你，帮我大忙了。。。回复 4# rubyish

作者: mcshell 时间: 2013-03-13 18:50

#!/usr/bin/perl
my (%hash,@dna,$text);
my $data = do {local $/;<DATA>};
my @line = split/^\s+/m,$data;
for (@line){
if(/([ATGC]+)[\s\S]+?([ATGC]+)/){
my ($len,$count) = length $2;
my @arr = split//,$2;
for (split//,$1){
$count++ if($_ eq shift(@arr));
}
my $m=$count/$len;
my $text = $m<0.6?1:$m<0.9?2:3;
open FH,">>$text.txt" or die "$!";
print FH $_;
close FH;
}
}
__DATA__
I001 195 TTGGT 205 225
E012 148 TTCGT 158 200
I001 1124 TTTCTT 1151 2978
E012 4 TAATGT 31 128
I002 17 AATTGACT 29 54
E023 5 AATTGACT 17 128
I003 16 TCGACTTAAT 35 65
E034 44 TGGACATAAT 63 200

复制代码

回复 1# zhagnqiang829

作者: sinian126 时间: 2013-03-14 19:04
提示: 作者被禁止或删除内容自动屏蔽

作者: mcshell 时间: 2013-03-14 21:13
回复 9# sinian126

所有字符串包括换行

欢迎光临 Chinaunix (http://bbs.chinaunix.net/)