跪求高手幫助解決perl問(wèn)題

弦斷有誰(shuí)聽(tīng)1053476508

回復(fù) 10# huang6894


    那就是真實(shí)文件，只是一部分。源文件我上傳了，你看看。在命令行里輸入Perl XX.pl  路徑\Orthologs-msu.txt Allseq.fa   這樣測(cè)試，對(duì)不對(duì)？

huang6894

回復(fù) 11# 弦斷有誰(shuí)聽(tīng)1053476508

不好意思，幫不到你，我這里測(cè)試沒(méi)有問(wèn)題~
   

b114213903

請(qǐng)?zhí)峁┩暾�的fa文件

弦斷有誰(shuí)聽(tīng)1053476508

Allseq.fa文件有200g，太大沒(méi)法上傳，只能發(fā)測(cè)試文件。大神，我有重新上傳了測(cè)試文件，編輯了帖子，麻煩您給看看，或許是過(guò)程步驟我沒(méi)有向您描述清楚。回復(fù) 12# huang6894


   

弦斷有誰(shuí)聽(tīng)1053476508

Allseq.fa文件有200g，太大沒(méi)法上傳，只能發(fā)測(cè)試文件，大神，我新上傳了測(cè)試文件，您再看看。回復(fù) 13# b114213903


   

弦斷有誰(shuí)聽(tīng)1053476508

還請(qǐng)大家不吝賜教!

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";

復(fù)制代碼

回復(fù) 16# 弦斷有誰(shuí)聽(tīng)1053476508
明天給你第三步的

   

chenhao392

可能的小老鄉(xiāng)你好.

既然你是洛陽(yáng)上的本科，我猜我們多半 算半個(gè)老鄉(xiāng)。我也有高中時(shí)代的哥們?cè)谀撬�大學(xué)讀本科�？丛谖覀兪抢相l(xiāng)的份上，看在你不停的積極討論的份上，我也就不潛水了。程序大概給你寫(xiě)完了，請(qǐng)查收。
我理解你的需求是：根據(jù)OrthoMCL計(jì)算出的同源基因名稱(chēng)，從不同的植物的基因序列 (200 Gb)中，拿到OrthoMCL提到的基因，并把一個(gè)同源基因組的存放在一個(gè)文件里。


但是！
每一次提問(wèn)前，請(qǐng)一定要想一想自己的問(wèn)題是否問(wèn)的夠有效率。這次你忽略了幾個(gè)比較關(guān)鍵的信息：
1. Step 1-3 似乎是你自己的programming想法。其實(shí)不妨把整個(gè)需求描述出來(lái)，也許別人有更好的辦法呢？ 比如我的code，用到了hash of hash.

2. 一個(gè)基因會(huì)有多個(gè)isoforms，你的ortholog 信息是gene的，而你提供的Allseq.fa 是根據(jù)isoform的，直接的ID當(dāng)然一個(gè)match都沒(méi)有。我的code假設(shè)你的Allseq.fa 里的isoform是加了 . (點(diǎn)) 或者 _(下劃線(xiàn)) 的。這cover了玉米，水稻和擬南芥，并不一定cover了所有的情況。所以，這個(gè)代碼你可能還要改。


3. Again, 一個(gè)gene有多個(gè)isoforms，要選哪一個(gè)呢？這要看你的研究課題了。我知道通常會(huì)選最長(zhǎng)的那一個(gè)。

4. 請(qǐng)以后盡量學(xué)習(xí)用linux，你用windows 下載的Allseq.fa 的換行符，需要多用一個(gè)chop，將來(lái)記不住自己的代碼干了什么的話(huà)，會(huì)出錯(cuò)。

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

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

復(fù)制代碼

在我看來(lái)，數(shù)字416 是說(shuō)這個(gè)ortholog在這些物種里一共有多少個(gè)copy，數(shù)字 2是說(shuō)有多少個(gè)物種里有這個(gè)ortholog，在上文中是rice(水稻)和 sorghum(高粱)。我的程序已經(jīng)考慮了這一點(diǎn)，會(huì)從類(lèi)似 LOC_Os04g04140 的位置讀取。

請(qǐng)用嚴(yán)謹(jǐn)認(rèn)真的態(tài)度來(lái)做科研和在版面上提出問(wèn)題。將來(lái)perl玩熟練的話(huà)，也請(qǐng)回到版面上幫助別人。

        此致，
你的半個(gè)老鄉(xiāng)： chenhao392


run

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

復(fù)制代碼

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. }
    

復(fù)制代碼

弦斷有誰(shuí)聽(tīng)1053476508

太感謝了！大神，你怎么知道我在洛陽(yáng)上的啊，我在河科大上的本科。你寫(xiě)的東西我認(rèn)真看了，很受啟發(fā)呢。不嫌棄的話(huà)交個(gè)朋友吧，我q1053476508回復(fù) 18# chenhao392


   

b114213903

回復(fù) 18# chenhao392


    200G的fa文件，這樣讀，能行么？