求助：取最长转录本脚本

yang7473453

现有一个fasta文件，里面是人类编码区序列数据，内容如下两图


ENSG是基因编号，ENST是转录本编号，图中可见同一个基因编号有两个转录本，当然实际上可能一个基因编号可能有多个转录本，现需要一个基因编号只留最长的一条转录本，如何提取？求perl脚本

aku1

能贴出部分实际内容吗？

yang7473453

已经贴出来了呀！图中可看到>后有两个编号一个是ENSG的一个是ENST的，ENSG编号会有相同的，当ENSG编号相同时我们只取序列最长的那一条回复 2# aku1


   

b114213903

回复 3# yang7473453


    你给的图片，怎么做脚本测试？

    给个示例文件，很困难么？

yang7473453

懂了，我贴一段上来，下面这些序列中有两个基因，ENSG编号相同的我需要取最长的那条回复 4# b114213903
>ENSG00000104671|ENST00000522141
ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
AGTGAAATCAGAGGAGATGTAACTATCGTTACCCAGATAATATCACTCCTGACACTGAAG
ATCCAGAACCAAAACCTATGA
>ENSG00000104671|ENST00000520829
ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
CTTATCATAAATGCTTACCCAGATAATATCACTCCTGACACTGAAGATCCAGAACCAAAA
CCTATGATCATTGGCACCAATAATGTGTTTGAAGTTGGCTGTTATTCCCAAGCCATGAAG
ATGGGAGATAATAATGTCATTGAATCAAAAGCATATGTAGGCAGAAATGTAATATTGACA
AGTGGCTGCATCATTGGGGCTTGTTGCAACCTAAATACATTTGAAGTCATCCCTGAGAAT
ACGGTGATCTATGGTGCAGACTGCCTTCGTCGGGTGCAGACTGAGCGACCGCAGGTACTA
GAACCTCTCTTTAAAAAGAGTTCTATCTGCTGA
>ENSG00000104671|ENST00000523666
ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
CTTATCATAAATGCATTCCCAAGCCATGAAGATGGGAGATAA
>ENSG00000205339|ENST00000533233
ATGGACCCCAACACCATTATCGAGGCCCTGCGGGGCACTATGGACCCAGCCCTGCGTGAG
GCCGCGGAGCGCCAGCTCAATGAAACGGAGTTTCGCTCTTGTTGCCCAGACTGGAGTGCA
ATGGCACAATCTCGGCTCGCCACAACCTCCGCCTCCCGGATTCGGGCAATTCTCCTGCCT
CAGCCTCTTGAGTAG
>ENSG00000205339|ENST00000527431
ATGATAACACAGTATTGGCCTGATCGAGAAACAGCACCAGGGGATATATCCCCTTATACT
ATTCCAGAAGAAGATCGCCATTGTATTCGAGAAAATATTGTAGAAGCCATTATCCATTCT
CCTGAGCTCATCAGGGTACAGCTTACTACATGCATTCATCACATCATCAAACATGATTAT
CCAAGCCGCTGGACTGCCATTGTGGACAAAATTGGCTTTTATCTTCAGTCCGATAACAGT
GCTTGTTGGCTAGGAATTCTTCTTTGCCTTTATCAGCTTGTGAAAAATTATGAGTATAAA
AAACCAGAGGAGCGGAGTCCATTGGTAGCAGCAATGCAGCATTTTCTGCCAGTTCTAAAG
GATCGTTTTATCCAGCTTCTTTCTGACCAGTC

   

b114213903

回复 5# yang7473453

1. #!/usr/bin/perl
   
2. use Bio::SeqIO;
   
3. use strict;
   
4. 
   
5. ##需要提供序列文件##
   
6. unless (@ARGV==1){
   
7.         die "Usage:\t$0  seq_file_name\n";
   
8. }
   
9. 
   
10. ##读取并比较序列长短##
    
11. my $Fas=Bio::SeqIO->new(-file=>"$ARGV[0]",-format=>"fasta") or die "Open seq file failed! $!\n";
    
12. my (%Length,%Seqs)=();
    
13. while(my $seq=$Fas->next_seq()){
    
14.         my $id=$seq->id();
    
15.         my ($ensg,$enst)=split (/\|/,$id);
    
16.         unless($Length{$ensg}){
    
17.                 $Length{$ensg}=$seq->length();
    
18.                 $Seqs{$ensg}=$seq;
    
19.         }elsif($Length{$ensg} < $seq->length()){
    
20.                 $Length{$ensg}=$seq->length();
    
21.                 $Seqs{$ensg}=$seq;
    
22.         }
    
23. }
    
24. $Fas->close();
    
25. 
    
26. ##输出比较结果##
    
27. my $Out=Bio::SeqIO->new(-file=>">Out.fasta",-format=>"fasta") or die "Write to Out seq failed! $!\n";
    
28. foreach my $key(sort keys %Seqs){
    
29.         print "Now is:\t$key\n";
    
30.         $Out->write_seq($Seqs{$key});
    
31. }
    
32. $Out->close();
    
33. print "Finised!\n";

复制代码

1. >ENSG00000104671|ENST00000522141
   
2. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
   
3. AGTGAAATCAGAGGAGATGTAACTATCGTTACCCAGATAATATCACTCCTGACACTGAAG
   
4. ATCCAGAACCAAAACCTATGA
   
5. >ENSG00000104671|ENST00000520829
   
6. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
   
7. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
   
8. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
   
9. CTTATCATAAATGCTTACCCAGATAATATCACTCCTGACACTGAAGATCCAGAACCAAAA
   
10. CCTATGATCATTGGCACCAATAATGTGTTTGAAGTTGGCTGTTATTCCCAAGCCATGAAG
    
11. ATGGGAGATAATAATGTCATTGAATCAAAAGCATATGTAGGCAGAAATGTAATATTGACA
    
12. AGTGGCTGCATCATTGGGGCTTGTTGCAACCTAAATACATTTGAAGTCATCCCTGAGAAT
    
13. ACGGTGATCTATGGTGCAGACTGCCTTCGTCGGGTGCAGACTGAGCGACCGCAGGTACTA
    
14. GAACCTCTCTTTAAAAAGAGTTCTATCTGCTGA
    
15. >ENSG00000104671|ENST00000523666
    
16. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
17. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
    
18. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
    
19. CTTATCATAAATGCATTCCCAAGCCATGAAGATGGGAGATAA
    
20. >ENSG00000205339|ENST00000533233
    
21. ATGGACCCCAACACCATTATCGAGGCCCTGCGGGGCACTATGGACCCAGCCCTGCGTGAG
    
22. GCCGCGGAGCGCCAGCTCAATGAAACGGAGTTTCGCTCTTGTTGCCCAGACTGGAGTGCA
    
23. ATGGCACAATCTCGGCTCGCCACAACCTCCGCCTCCCGGATTCGGGCAATTCTCCTGCCT
    
24. CAGCCTCTTGAGTAG
    
25. >ENSG00000205339|ENST00000527431
    
26. ATGATAACACAGTATTGGCCTGATCGAGAAACAGCACCAGGGGATATATCCCCTTATACT
    
27. ATTCCAGAAGAAGATCGCCATTGTATTCGAGAAAATATTGTAGAAGCCATTATCCATTCT
    
28. CCTGAGCTCATCAGGGTACAGCTTACTACATGCATTCATCACATCATCAAACATGATTAT
    
29. CCAAGCCGCTGGACTGCCATTGTGGACAAAATTGGCTTTTATCTTCAGTCCGATAACAGT
    
30. GCTTGTTGGCTAGGAATTCTTCTTTGCCTTTATCAGCTTGTGAAAAATTATGAGTATAAA
    
31. AAACCAGAGGAGCGGAGTCCATTGGTAGCAGCAATGCAGCATTTTCTGCCAGTTCTAAAG
    
32. GATCGTTTTATCCAGCTTCTTTCTGACCAGTC

复制代码

1. >ENSG00000104671|ENST00000520829
   
2. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
   
3. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
   
4. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
   
5. CTTATCATAAATGCTTACCCAGATAATATCACTCCTGACACTGAAGATCCAGAACCAAAA
   
6. CCTATGATCATTGGCACCAATAATGTGTTTGAAGTTGGCTGTTATTCCCAAGCCATGAAG
   
7. ATGGGAGATAATAATGTCATTGAATCAAAAGCATATGTAGGCAGAAATGTAATATTGACA
   
8. AGTGGCTGCATCATTGGGGCTTGTTGCAACCTAAATACATTTGAAGTCATCCCTGAGAAT
   
9. ACGGTGATCTATGGTGCAGACTGCCTTCGTCGGGTGCAGACTGAGCGACCGCAGGTACTA
   
10. GAACCTCTCTTTAAAAAGAGTTCTATCTGCTGA
    
11. >ENSG00000205339|ENST00000527431
    
12. ATGATAACACAGTATTGGCCTGATCGAGAAACAGCACCAGGGGATATATCCCCTTATACT
    
13. ATTCCAGAAGAAGATCGCCATTGTATTCGAGAAAATATTGTAGAAGCCATTATCCATTCT
    
14. CCTGAGCTCATCAGGGTACAGCTTACTACATGCATTCATCACATCATCAAACATGATTAT
    
15. CCAAGCCGCTGGACTGCCATTGTGGACAAAATTGGCTTTTATCTTCAGTCCGATAACAGT
    
16. GCTTGTTGGCTAGGAATTCTTCTTTGCCTTTATCAGCTTGTGAAAAATTATGAGTATAAA
    
17. AAACCAGAGGAGCGGAGTCCATTGGTAGCAGCAATGCAGCATTTTCTGCCAGTTCTAAAG
    
18. GATCGTTTTATCCAGCTTCTTTCTGACCAGTC
    

复制代码

yang7473453

非常感谢，写的很好，很好用。回复 6# b114213903


   

huang6894

1. #!/usr/bin/perl -w
   
2. my(%h1, %h2);
   
3. $/="\>";
   
4. <DATA>;
   
5. while(<DATA>){
   
6.         chomp;
   
7.         my($g_t, @fastq) = (split /\n/);
   
8.         my($gene, $tran) = (split /\|/,$g_t);
   
9.         push @{$h2{$tran}},@fastq;
   
10.         $h1{$gene} ||= $tran ;
    
11.         $h1{$gene} = $tran  if length(join("",@{$h2{$h1{$gene}}})) < length(join("",@fastq));       
    
12. }
    
13. $/="\n";
    
14. foreach my $gene(keys %h1){
    
15.         print ">$gene|$h1{$gene}\n";
    
16.         print join("\n",@{$h2{$h1{$gene}}});
    
17.         print $/;
    
18. }
    
19. __DATA__
    
20. >ENSG00000104671|ENST00000522141
    
21. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
22. AGTGAAATCAGAGGAGATGTAACTATCGTTACCCAGATAATATCACTCCTGACACTGAAG
    
23. ATCCAGAACCAAAACCTATGA
    
24. >ENSG00000104671|ENST00000520829
    
25. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
26. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
    
27. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
    
28. CTTATCATAAATGCTTACCCAGATAATATCACTCCTGACACTGAAGATCCAGAACCAAAA
    
29. CCTATGATCATTGGCACCAATAATGTGTTTGAAGTTGGCTGTTATTCCCAAGCCATGAAG
    
30. ATGGGAGATAATAATGTCATTGAATCAAAAGCATATGTAGGCAGAAATGTAATATTGACA
    
31. AGTGGCTGCATCATTGGGGCTTGTTGCAACCTAAATACATTTGAAGTCATCCCTGAGAAT
    
32. ACGGTGATCTATGGTGCAGACTGCCTTCGTCGGGTGCAGACTGAGCGACCGCAGGTACTA
    
33. GAACCTCTCTTTAAAAAGAGTTCTATCTGCTGA
    
34. >ENSG00000104671|ENST00000523666
    
35. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
36. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
    
37. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
    
38. CTTATCATAAATGCATTCCCAAGCCATGAAGATGGGAGATAA
    
39. >ENSG00000205339|ENST00000533233
    
40. ATGGACCCCAACACCATTATCGAGGCCCTGCGGGGCACTATGGACCCAGCCCTGCGTGAG
    
41. GCCGCGGAGCGCCAGCTCAATGAAACGGAGTTTCGCTCTTGTTGCCCAGACTGGAGTGCA
    
42. ATGGCACAATCTCGGCTCGCCACAACCTCCGCCTCCCGGATTCGGGCAATTCTCCTGCCT
    
43. CAGCCTCTTGAGTAG
    
44. >ENSG00000205339|ENST00000527431
    
45. ATGATAACACAGTATTGGCCTGATCGAGAAACAGCACCAGGGGATATATCCCCTTATACT
    
46. ATTCCAGAAGAAGATCGCCATTGTATTCGAGAAAATATTGTAGAAGCCATTATCCATTCT
    
47. CCTGAGCTCATCAGGGTACAGCTTACTACATGCATTCATCACATCATCAAACATGATTAT
    
48. CCAAGCCGCTGGACTGCCATTGTGGACAAAATTGGCTTTTATCTTCAGTCCGATAACAGT
    
49. GCTTGTTGGCTAGGAATTCTTCTTTGCCTTTATCAGCTTGTGAAAAATTATGAGTATAAA
    
50. AAACCAGAGGAGCGGAGTCCATTGGTAGCAGCAATGCAGCATTTTCTGCCAGTTCTAAAG
    
51. GATCGTTTTATCCAGCTTCTTTCTGACCAGTC

复制代码

251744647

#awk 'NR==1{substr($0,2)}1' urfile |awk 'BEGIN{FS="|";RS="\n>"}{if(a[$1]<length($2)){a[$1]=length($2);b[$1]=$2}}END{for(i in a)print">" i"|"b[i]}' |sed '/^$/d'
>ENSG00000205339|ENST00000527431
ATGATAACACAGTATTGGCCTGATCGAGAAACAGCACCAGGGGATATATCCCCTTATACT
ATTCCAGAAGAAGATCGCCATTGTATTCGAGAAAATATTGTAGAAGCCATTATCCATTCT
CCTGAGCTCATCAGGGTACAGCTTACTACATGCATTCATCACATCATCAAACATGATTAT
CCAAGCCGCTGGACTGCCATTGTGGACAAAATTGGCTTTTATCTTCAGTCCGATAACAGT
GCTTGTTGGCTAGGAATTCTTCTTTGCCTTTATCAGCTTGTGAAAAATTATGAGTATAAA
AAACCAGAGGAGCGGAGTCCATTGGTAGCAGCAATGCAGCATTTTCTGCCAGTTCTAAAG
GATCGTTTTATCCAGCTTCTTTCTGACCAGTC
>ENSG00000104671|ENST00000520829
ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
CTTATCATAAATGCTTACCCAGATAATATCACTCCTGACACTGAAGATCCAGAACCAAAA
CCTATGATCATTGGCACCAATAATGTGTTTGAAGTTGGCTGTTATTCCCAAGCCATGAAG
ATGGGAGATAATAATGTCATTGAATCAAAAGCATATGTAGGCAGAAATGTAATATTGACA
AGTGGCTGCATCATTGGGGCTTGTTGCAACCTAAATACATTTGAAGTCATCCCTGAGAAT
ACGGTGATCTATGGTGCAGACTGCCTTCGTCGGGTGCAGACTGAGCGACCGCAGGTACTA
GAACCTCTCTTTAAAAAGAGTTCTATCTGCTGA

清泉一边

一段时间没写了，感觉好多都忘了啊，试试吧~

1. #!/usr/bin/perl  -w
   
2. use strict;
   
3. my $ref = {};
   
4. $/ = ">";<DATA>;$/= "\n";
   
5. while (<DATA>){
   
6. if ($_ =~ /(\S+)\|(\S+)/){
   
7. my $gene = $1;
   
8. my $gi = $2;
   
9. $/ = ">";
   
10. my $seq = <DATA>;
    
11. $seq =~ s/>$//;
    
12. $ref->{$gene}->{$seq} = $gi;
    
13. }
    
14. $/ = "\n";
    
15. }
    
16. foreach my $gene (sort keys %$ref){
    
17. my $max = 0;
    
18. my $max_length = 0;
    
19.   foreach my $seq (sort keys $ref->{$gene}){
    
20.    if ($max_length < length $seq){
    
21.     $max = $seq;
    
22.     $max_length =length $seq;
    
23.      }
    
24.   }
    
25.   print ">$gene|$ref->{$gene}->{$max}\n$max\n" if exists $ref->{$gene}->{$max};
    
26. }
    
27. 
    
28. __DATA__
    
29. >ENSG00000104671|ENST00000522141
    
30. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
31. AGTGAAATCAGAGGAGATGTAACTATCGTTACCCAGATAATATCACTCCTGACACTGAAG
    
32. ATCCAGAACCAAAACCTATGA
    
33. >ENSG00000104671|ENST00000520829
    
34. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
35. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
    
36. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
    
37. CTTATCATAAATGCTTACCCAGATAATATCACTCCTGACACTGAAGATCCAGAACCAAAA
    
38. CCTATGATCATTGGCACCAATAATGTGTTTGAAGTTGGCTGTTATTCCCAAGCCATGAAG
    
39. ATGGGAGATAATAATGTCATTGAATCAAAAGCATATGTAGGCAGAAATGTAATATTGACA
    
40. AGTGGCTGCATCATTGGGGCTTGTTGCAACCTAAATACATTTGAAGTCATCCCTGAGAAT
    
41. ACGGTGATCTATGGTGCAGACTGCCTTCGTCGGGTGCAGACTGAGCGACCGCAGGTACTA
    
42. GAACCTCTCTTTAAAAAGAGTTCTATCTGCTGA
    
43. >ENSG00000104671|ENST00000523666
    
44. ATGGCGGAGAAGACTCAAAAGAGTGTGAAGATTGCTCCTGGAGCAGTTGTATGTGTAGAA
    
45. AGTGAAATCAGAGGAGATGTAACTATCGGACCTCGGACAGTGATCCACCCTAAAGCAAGA
    
46. ATTATTGCGGAAGCCGGGCCAATAGTGATTGGCGAAGGGAACCTAATAGAAGAACAGGCC
    
47. CTTATCATAAATGCATTCCCAAGCCATGAAGATGGGAGATAA
    
48. >ENSG00000205339|ENST00000533233
    
49. ATGGACCCCAACACCATTATCGAGGCCCTGCGGGGCACTATGGACCCAGCCCTGCGTGAG
    
50. GCCGCGGAGCGCCAGCTCAATGAAACGGAGTTTCGCTCTTGTTGCCCAGACTGGAGTGCA
    
51. ATGGCACAATCTCGGCTCGCCACAACCTCCGCCTCCCGGATTCGGGCAATTCTCCTGCCT
    
52. CAGCCTCTTGAGTAG
    
53. >ENSG00000205339|ENST00000527431
    
54. ATGATAACACAGTATTGGCCTGATCGAGAAACAGCACCAGGGGATATATCCCCTTATACT
    
55. ATTCCAGAAGAAGATCGCCATTGTATTCGAGAAAATATTGTAGAAGCCATTATCCATTCT
    
56. CCTGAGCTCATCAGGGTACAGCTTACTACATGCATTCATCACATCATCAAACATGATTAT
    
57. CCAAGCCGCTGGACTGCCATTGTGGACAAAATTGGCTTTTATCTTCAGTCCGATAACAGT
    
58. GCTTGTTGGCTAGGAATTCTTCTTTGCCTTTATCAGCTTGTGAAAAATTATGAGTATAAA
    
59. AAACCAGAGGAGCGGAGTCCATTGGTAGCAGCAATGCAGCATTTTCTGCCAGTTCTAAAG
    
60. GATCGTTTTATCCAGCTTCTTTCTGACCAGTC

复制代码

回复 1# yang7473453