Chinaunix

标题: perl [打印本页]

---

作者: super_two    时间: 2015-08-07 21:07
标题: perl
我的文件 都是这种FASTA格式：
qwer.1         AFRTYUIFGHJKLVBNM
qwer.2         BVXNVFGSFYEBCSHB
qwer.3         HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC
rtyuip00.1    AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
hhjkl.1         JDNVKHFBAKHFAKFAFJNSADFJAS
hhjk.2         HFSDHNFKANFAKFIJI
hhjk.3         JNFAJNFALSDFLMAD
hhjk.4         KJGSEGJOAKFFDSMFAPOKEF
...
现在我想把小数点前相同的放到一起，点后的忽略。然后比较小数点前相同的ID后的序列长度，将最长的输出来；若有多条一样长则随机输出其中一条，请大神用perl写个程序。其实就是处理可变剪切，不知道我是否描述清楚了。

---

作者: super_two    时间: 2015-08-07 21:10
标题: perl 匹配
我的文件 都是这种FASTA格式：
qwer.1         AFRTYUIFGHJKLVBNM
qwer.2         BVXNVFGSFYEBCSHB
qwer.3         HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC
rtyuip00.1    AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
hhjkl.1         JDNVKHFBAKHFAKFAFJNSADFJAS
hhjk.2         HFSDHNFKANFAKFIJI
hhjk.3         JNFAJNFALSDFLMAD
hhjk.4         KJGSEGJOAKFFDSMFAPOKEF
...
现在我想把小数点前相同的放到一起，点后的忽略。然后比较小数点前相同的ID后的序列长度，将最长的输出来；若有多条一样长则随机输出其中一条，请大神用perl写个程序。其实就是处理可变剪切，不知道我是否描述清楚了。

---

作者: MMMIX    时间: 2015-08-07 21:55
@b114213903

---

作者: zxy877298415    时间: 2015-08-07 22:25

1. awk -F '[. ]+' '{a[$1]=a[$1]?a[$1]" "$3:$3}END{for(i in a){printf i" ";for(j=1;j<=split(a[i],b," ");j++) {if(length(b[j])>int(t))t=b[j]}print t;t=""}}' file

复制代码

回复 1# super_two


   

---

作者: super_two    时间: 2015-08-07 22:31
我是新手，麻烦您给解释一下，谢谢~回复 2# zxy877298415


   

---

作者: tc1989tc    时间: 2015-08-07 22:32
字典加数组 然后sort数组

---

作者: MMMIX    时间: 2015-08-07 22:45
为啥不自己先写一个呢? 新手就应该多练习嘛

---

作者: super_two    时间: 2015-08-07 23:12
我写了一个，但是不能运行。我的想法就是用Bio::SeqIO模块，提取出原始文件的ID，和序列，并以内length计算出每个序列的长度。
将提出的ID按照小数点前相同的ID及其序列放进一个数组，然后分别在这个数组中找出最长的序列输出来，但是要将原始ID、序列以及长度都分别对应，我不知道如何才能对应好。不知道我这个思路对不对




回复 7# MMMIX


   

---

作者: MMMIX    时间: 2015-08-08 00:29
回复 8# super_two


    不能运行也罢, 贴上来看看.

---

作者: 聆雨淋夜    时间: 2015-08-08 00:31


想得复杂点了

1. awk -F '[. ]+' '$1!=a{ml=length($3);a=$1;b[$1]=$3}$1==a{if(length($3)>ml){ml=length($3);b[$1]=$3}}END{for(i in b)print i"  "b[i]}' file

复制代码

hhjk  KJGSEGJOAKFFDSMFAPOKEF
hhjkl  JDNVKHFBAKHFAKFAFJNSADFJAS
rtyuip00  AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
qwer  HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC

---

作者: super_two    时间: 2015-08-08 10:40
这是我写的，结果不对，我并没有办法将我说的思路用perl写出来，麻烦您看看
#!/usr/bin/perl   -w
use warnings;
use strict;
use Bio::SeqIO;

my $data=$ARGV[0];

open (W,"<$data") || die $!;
    my %gene;
    while(<W>){
    chomp;
    $gene{$_}++;
}
close W;

open (M,"> $data.longest") || die $!;

my $all=Bio::SeqIO->new(-file=>"$data",-format=>'fasta');
while(my $seq=$all->next_seq){
    my $id=$seq->id;
    my $seq=$seq->seq;
     $gene{$seq}=length($seq);
push(my @length,$gene{$seq});
my @ids=split(/\|/,$id);
    my @newid=split(/\./,$id);
    foreach $a (@newid){
        if ($a){
            my @length_1 = sort{$a<=>$b}@length;
            print ">$id\n";
            print "$seq\n";
        }
        else{
            print ">$id\t";
            print "$seq\n";


        }
    }
}
close M;回复 9# MMMIX


   

---

作者: super_two    时间: 2015-08-08 11:12
这是我写的，结果不对，我并没有办法将我说的思路用perl写出来，麻烦您看看
#!/usr/bin/perl   -w
use warnings;
use strict;
use Bio::SeqIO;

my $data=$ARGV[0];

open (W,"<$data") || die $!;
    my %gene;
    while(<W>){
    chomp;
    $gene{$_}++;
}
close W;

open (M,"> $data.longest") || die $!;

my $all=Bio::SeqIO->new(-file=>"$data",-format=>'fasta');
while(my $seq=$all->next_seq){
    my $id=$seq->id;
    my $seq=$seq->seq;
     $gene{$seq}=length($seq);
push(my @length,$gene{$seq});
my @ids=split(/\|/,$id);
    my @newid=split(/\./,$id);
    foreach $a (@newid){
        if ($a){
            my @length_1 = sort{$a<=>$b}@length;
            print ">$id\n";
            print "$seq\n";
        }
        else{
            print ">$id\t";
            print "$seq\n";


        }
    }
}
close M

---

作者: j_cle    时间: 2015-08-08 11:53
回复 10# 聆雨淋夜

1. awk -F '[. ]+' '$1!=a{ml=length($3);a=$1;b[$1]=$3}$1==a{if(length($3)>ml){ml=length($3);b[$1]=$3}}END{for(i in b)print i"  "b[i]}' file

复制代码

请教下 $1!=a{}$1==a{}  这个结构是什么意思？以前还没有见过这种写法，
我自己猜想是
如果在匹配的过程中   $1如果相同的话执行==后面的表达式{}，如果$1不相同的话执行！=后面的表达式
是这样理解的么？
这种用法是叫什么呢？哪里能找到相关的资料，想具体学习下，谢谢！

---

作者: MMMIX    时间: 2015-08-08 14:31
回复 11# super_two


    其实, 你只需要逐个读取序列, 对于每个 id, 只保留最长的那个序列即可.

---

作者: 聆雨淋夜    时间: 2015-08-08 15:23
回复 13# j_cle

a是上一行的$1，$1!=a 表示当前行的$1和上一行的$1不一样了，这样就进行ml和b[$1]的初始化


   

---

作者: 聆雨淋夜    时间: 2015-08-08 15:25
回复 13# j_cle
其实不需要这么区分的，见CU的大神写过类似的，只是我想不起来怎么写的了。


   

---

作者: super_two    时间: 2015-08-08 17:49
您可以写一个脚本让我看一下整个语法过程吗？见过的脚本太少，无法把握整体的语法结构 回复 15# 聆雨淋夜


   

---

作者: super_two    时间: 2015-08-08 18:54
您可以写一个脚本让我看一下整个语法过程吗？见过的脚本太少，无法把握整体的语法结构回复 14# MMMIX


   

---

作者: MMMIX    时间: 2015-08-08 19:19


回复 18# super_two


    真心建议你找些 Perl 编程的书好好看看, 不过估计你是听不进去的.

无论如何, 下述脚本演示了解决你的问题的一种方法. 另外,如果在你的输入文件中,id相同的序列顺序出现,那么也可以不用hash.

1. 
   
2. use strict;
   
3. use warnings;
   
4. 
   
5. use v5.14;
   
6. 
   
7. my %record;
   
8. while (<DATA>) {
   
9.     my @f = split;
   
10.     $f[0] =~ s/(?:\.\d*)?$//;
    
11.     if (! $record{$f[0]} || $record{$f[0]}->[0] < length($f[1])) {
    
12.         $record{$f[0]} = [length $f[1], $f[1]];
    
13.     }
    
14. }
    
15. 
    
16. say "output version 1:";
    
17. while (my ($id, $seq) = each %record) {
    
18.     say "$id\t", $seq->[1];
    
19. }
    
20. 
    
21. say "output version 2:";
    
22. for my $id (sort keys %record) {
    
23.     printf "%-12s %s\n", $id, $record{$id}->[1];
    
24. }
    
25. 
    
26. __DATA__
    
27. qwer.1         AFRTYUIFGHJKLVBNM
    
28. qwer.2         BVXNVFGSFYEBCSHB
    
29. qwer.3         HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC
    
30. rtyuip00.1    AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
    
31. hhjkl.1         JDNVKHFBAKHFAKFAFJNSADFJAS
    
32. hhjk.2         HFSDHNFKANFAKFIJI
    
33. hhjk.3         JNFAJNFALSDFLMAD
    
34. hhjk.4         KJGSEGJOAKFFDSMFAPOKEF
    

复制代码

---

作者: 聆雨淋夜    时间: 2015-08-08 19:37
回复 19# MMMIX
perl可读性好差


   

---

作者: MMMIX    时间: 2015-08-08 21:19
回复 20# 聆雨淋夜


    对于我不熟悉的语言来说, 我就没有发现可读性好的.

---

作者: b114213903    时间: 2015-08-12 08:43


回复 3# MMMIX

1. #!/usr/bin/perl
   
2. use strict;
   
3. use Bio::SeqIO;
   
4. 
   
5. my $fasta=shift @ARGV;
   
6. (my $Out=$fasta)=~s/(\.[^\.]+$)/_out$1/;
   
7. 
   
8. my $IN=Bio::SeqIO->new(-file=>"$fasta",-format=>'fasta');
   
9. my $Lenths={};
   
10. my $Seqs={};
    
11. while(my $seq=$IN->next_seq){
    
12.         my $id=$seq->id;
    
13.         my $length=$seq->length;
    
14.         (my $new_id=$id)=~s/^([^\.]+)\..+$/$1/;
    
15.        
    
16.         if($Lenths->{$new_id}){
    
17.                 if($Lenths->{$new_id}<$length){
    
18.                         $Lenths->{$new_id}=$length;
    
19.                         $Seqs->{$new_id}=$seq;
    
20.                 }
    
21.         }else{
    
22.                 $Lenths->{$new_id}=$length;
    
23.                 $Seqs->{$new_id}=$seq;
    
24.         }
    
25. }
    
26. $IN->close();
    
27. 
    
28. my $OUT=Bio::SeqIO->new(-file=>">$Out",-format=>'fasta');
    
29. foreach my $id(sort keys %{$Seqs}){
    
30.         print "$id\n";
    
31.         $OUT->write_seq($Seqs->{$id});
    
32. }
    
33. $OUT->close();

复制代码

1. >qwer.1
   
2. AFRTYUIFGHJKLVBNM
   
3. >qwer.2
   
4. BVXNVFGSFYEBCSHB
   
5. >qwer.3
   
6. HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC
   
7. >rtyuip00.1
   
8. AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
   
9. >hhjkl.1
   
10. JDNVKHFBAKHFAKFAFJNSADFJAS
    
11. >hhjk.2
    
12. HFSDHNFKANFAKFIJI
    
13. >hhjk.3
    
14. JNFAJNFALSDFLMAD
    
15. >hhjk.4
    
16. KJGSEGJOAKFFDSMFAPOKEF

复制代码

1. >hhjk.4
   
2. KJGSEGJOAKFFDSMFAPOKEF
   
3. >hhjkl.1
   
4. JDNVKHFBAKHFAKFAFJNSADFJAS
   
5. >qwer.3
   
6. HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC
   
7. >rtyuip00.1
   
8. AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
   

复制代码

另外，FASTA格式的文档，有规定格式：一、ID行：“>”（大于号）开头，后面根ID名称及各种注释信息，不可换行。二、数据行：单纯的核酸或氨基酸序列，可以单行也可多行。

---

作者: b114213903    时间: 2015-08-12 08:47
回复 8# super_two


    试试我用Bio::SeqIO模块的脚本！

---

作者: super_two    时间: 2015-08-12 18:10

1. qwer
   
2. AFRTYUIFGHJKLVBNM
   
3. qwer
   
4. AFRTYUIFGHJKLVBNM
   
5. qwer
   
6. HDFKSHFGSERYFIEURHFSUFDSHVBSJEUABFUHABFCAHFBC
   
7. qwer
   
8. AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
   
9. rtyuip00
   
10. AFHBVSFHUACFKUSHDBAKFHAKUFHSADKFUA
    
11. hhjkl
    
12. JDNVKHFBAKHFAKFAFJNSADFJAS
    
13. qwer
    
14. JDNVKHFBAKHFAKFAFJNSADFJAS
    
15. rtyuip00
    
16. JDNVKHFBAKHFAKFAFJNSADFJAS
    
17. hhjk
    
18. HFSDHNFKANFAKFIJI
    
19. hhjkl
    
20. HFSDHNFKANFAKFIJI
    
21. qwer
    
22. HFSDHNFKANFAKFIJI
    
23. rtyuip00
    
24. HFSDHNFKANFAKFIJI
    
25. hhjk
    
26. HFSDHNFKANFAKFIJI
    
27. hhjkl
    
28. HFSDHNFKANFAKFIJI
    
29. qwer
    
30. HFSDHNFKANFAKFIJI
    
31. rtyuip00
    
32. HFSDHNFKANFAKFIJI
    
33. hhjk
    
34. KJGSEGJOAKFFDSMFAPOKEF
    
35. hhjkl
    
36. KJGSEGJOAKFFDSMFAPOKEF
    
37. qwer
    
38. KJGSEGJOAKFFDSMFAPOKEF
    
39. rtyuip00
    
40. KJGSEGJOAKFFDSMFAPOKEF
    
41. 运行之后ID的小数点以及他后面的数字没有了，结果如上

复制代码

回复 23# b114213903


   

---

作者: b114213903    时间: 2015-08-12 19:55
回复 24# super_two


    我的结果里面都有啊

---

作者: MMMIX    时间: 2015-08-12 20:01
回复 24# super_two


    你是修改过人家的脚本吧?

---

作者: super_two    时间: 2015-08-13 11:27

1. hhjk
   
2. hhjkl
   
3. qwer
   
4. rtyuip00
   
5. 
   
6. 
   
7. 昨天我的编辑器保存出错，我完全复制粘贴上去的结果

复制代码

回复 26# MMMIX


   

---

作者: MMMIX    时间: 2015-08-13 11:35
回复 27# super_two


    你读懂那个脚本了么? 懂了的话就不会有疑问了.

---

作者: b114213903    时间: 2015-08-13 12:55
回复 26# MMMIX


    应该是有修改过，因为我用来保存数据的哈希引用中，都是由新ID指向的原始读入的序列，不可能会造成数据差异。

---

作者: MMMIX    时间: 2015-08-13 13:48
回复 29# b114213903


    我觉得他是光看了脚本的输出, 但并没有检查脚本写入的文件.

---

作者: super_two    时间: 2015-08-14 15:48
希望你们没被我蠢到，就算蠢到了，也别不理我哈我就是没看懂

1. (my $Out=$fasta)=~s/(\.[^\.]+$)/_out$1/;尤其是这句。

复制代码

也没看输入文件。@ b114213903@ MMMIX

---

作者: jason680    时间: 2015-08-14 15:56


回复 31# super_two

(my $Out=$fasta)=~s/(\.[^\.]+$)/_out$1/;  #尤其是这句。

as same as below:

my $Out = $fasta;   
$Out =~ s/(\.[^\.]+$)/_out$1/;   # 只 修改/替换 $Out
   

---

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

Powered by Discuz! X3.2