跪求高手帮助解决perl问题

弦断有谁听1053476508

回复 10# huang6894


    那就是真实文件，只是一部分。源文件我上传了，你看看。在命令行里输入Perl XX.pl  路径\Orthologs-msu.txt Allseq.fa   这样测试，对不对？

huang6894

回复 11# 弦断有谁听1053476508

不好意思，帮不到你，我这里测试没有问题~
   

b114213903

请提供完整的fa文件

弦断有谁听1053476508

Allseq.fa文件有200g，太大没法上传，只能发测试文件。大神，我有重新上传了测试文件，编辑了帖子，麻烦您给看看，或许是过程步骤我没有向您描述清楚。回复 12# huang6894


   

弦断有谁听1053476508

Allseq.fa文件有200g，太大没法上传，只能发测试文件，大神，我新上传了测试文件，您再看看。回复 13# b114213903


   

弦断有谁听1053476508

还请大家不吝赐教!

b114213903

1. #!/usr/bin/perl
   
2. use strict;
   
3. 
   
4. my $file=shift @ARGV || "orthol-msu.txt";
   
5. 
   
6. open (IN,"<$file") or die "Open $file failed!\n";
   
7. my %hash=();
   
8. while (my $line=<IN>){
   
9.         my @line=split /\t/,$line;
   
10.         my $value=$line[4];
    
11.         $line[4]=~s/\s/\n/g;
    
12.         print "$line[0]\n";
    
13.        
    
14.         ##第一步##
    
15.         open (OUT,">$line[0].txt");
    
16.         foreach my $item(@line){
    
17.                 print OUT "$item\n";
    
18.         }
    
19.         close OUT;
    
20.        
    
21.         ##第二步##
    
22.         $hash{$line[0]}=$value;
    
23. }
    
24. close IN;
    
25. foreach my $id(keys %hash){
    
26.         print "ID:\t$id\nValue:\t$hash{$id}\n";
    
27. }
    
28. print "Finished!\n";

复制代码

回复 16# 弦断有谁听1053476508
明天给你第三步的

   

chenhao392

可能的小老乡你好.

既然你是洛阳上的本科，我猜我们多半 算半个老乡。我也有高中时代的哥们在那所大学读本科。看在我们是老乡的份上，看在你不停的积极讨论的份上，我也就不潜水了。程序大概给你写完了，请查收。
我理解你的需求是：根据OrthoMCL计算出的同源基因名称，从不同的植物的基因序列 (200 Gb)中，拿到OrthoMCL提到的基因，并把一个同源基因组的存放在一个文件里。


但是！
每一次提问前，请一定要想一想自己的问题是否问的够有效率。这次你忽略了几个比较关键的信息：
1. Step 1-3 似乎是你自己的programming想法。其实不妨把整个需求描述出来，也许别人有更好的办法呢？ 比如我的code，用到了hash of hash.

2. 一个基因会有多个isoforms，你的ortholog 信息是gene的，而你提供的Allseq.fa 是根据isoform的，直接的ID当然一个match都没有。我的code假设你的Allseq.fa 里的isoform是加了 . (点) 或者 _(下划线) 的。这cover了玉米，水稻和拟南芥，并不一定cover了所有的情况。所以，这个代码你可能还要改。


3. Again, 一个gene有多个isoforms，要选哪一个呢？这要看你的研究课题了。我知道通常会选最长的那一个。

4. 请以后尽量学习用linux，你用windows 下载的Allseq.fa 的换行符，需要多用一个chop，将来记不住自己的代码干了什么的话，会出错。

5. 最重要的是，你真的理解你提供的文件么？
       Quote

1. 其中要求生成的ID文件如下：
   
2. APK_ORTHOMCL0
   
3. 416
   
4. 2
   
5. rice sorghum
   
6. LOC_Os04g04140
   
7. LOC_Os05g29160
   
8. LOC_Os10g02640

复制代码

在我看来，数字416 是说这个ortholog在这些物种里一共有多少个copy，数字 2是说有多少个物种里有这个ortholog，在上文中是rice(水稻)和 sorghum(高粱)。我的程序已经考虑了这一点，会从类似 LOC_Os04g04140 的位置读取。

请用严谨认真的态度来做科研和在版面上提出问题。将来perl玩熟练的话，也请回到版面上帮助别人。

        此致，
你的半个老乡： chenhao392


run

1. perl orthoMCL_picker.pl orthol-msu.txt Allseq.fa

复制代码

code

1. #!/usr/bin/perl
   
2. 
   
3. use strict;
   
4. use warnings;
   
5. 
   
6. my %id;
   
7. my %seq;
   
8. 
   
9. #load orthoMCL result from MSU
   
10. open FILE, "<$ARGV[0]" or die "$!\n";
    
11. while(<FILE>){
    
12.     chomp;
    
13.     #skip header line
    
14.     if($. == 1){
    
15.         next;
    
16.     }
    
17. 
    
18.     #get ortholog groups
    
19.     my ($orthoID,undef,$speciesCount,@line)=split(/\s+/,$_);
    
20.     for(my $i=$speciesCount;$i<scalar(@line);$i++){
    
21.         $id{$orthoID}{$line[$i]}="";
    
22.         $seq{$line[$i]}="";
    
23.     }
    
24. 
    
25. }
    
26. close FILE;
    
27. 
    
28. 
    
29. #scanning through the seq
    
30. my $switch=0;
    
31. my $id="";
    
32. open FILE, "<$ARGV[1]" or die "$!\n";
    
33. while(<FILE>){
    
34.     chomp;
    
35.     chop;
    
36.     if($_ =~ />(.*?)\s.*/){
    
37.         $switch =0;
    
38.         $id=$1;
    
39.         $id=~s/[\._].*//;
    
40.         if(defined $seq{$id}){
    
41.             $switch = 1;
    
42.         }
    
43.     }
    
44.     elsif($_ =~ /^[ATCGatcg]/ && $switch == 1){
    
45.         $seq{$id}.=$_;
    
46.     }
    
47. 
    
48. }
    
49. close FILE;
    
50. 
    
51. foreach my $orthoID(keys %id){
    
52.     open OUT, ">$orthoID.fa" or die "$!\n";
    
53.     foreach my $id(keys %{$id{$orthoID}}){
    
54.         print OUT ">$id\n$seq{$id}\n";
    
55.     }
    
56.     close OUT;
    
57. }
    

复制代码

弦断有谁听1053476508

太感谢了！大神，你怎么知道我在洛阳上的啊，我在河科大上的本科。你写的东西我认真看了，很受启发呢。不嫌弃的话交个朋友吧，我q1053476508回复 18# chenhao392


   

b114213903

回复 18# chenhao392


    200G的fa文件，这样读，能行么？