perl練習——FASTA格式文件中序列GC含量計算&perl數組排序如何獲得下標或者鍵

一、關於程序：

FUN：計算FASTA文件中每條序列中G和C的含量百分比，輸出最大值及其id

INPUT：FASTA格式文件

>seq1
CGCCGAGCGCTTGACCTCCAGCAAGACGCCGTCTGGCACATGCAACGAGCTGTAGCAGAC
>seq2
ATGCCTAGAACGTTCGAGACTTCTCGGGTGCGGTAGAATTAGCCATTCGACCGACTTCCA
GCATCTGCGAGCCGCCTGTTGATTGCATCCGCCGGGGACGCAACAAGGCAAGGCCCTAAC

OUTPUT：最高含量的序列id及其含量（這是上面的結果）

seq1
63.333333%

 

二、編程思想及代碼

 當是注釋行時（>……），獲得序列 ID ,並跳過該次循環；當讀到非注釋行即序列行時，記錄該行“G和C的含量”以及“序列的總含量”，這都可以利用perl上下文實現。（但是在這里有一些疑惑——當把14行@num換成$num會出現計算錯誤，知道的朋友歡迎留言）

use strict;
my %GC_content; # id=>GC_content
my %sequences; # id=>sequence
my ($id, $sum); # id, 每個序列的字符個數
my @num; # 中間變量，用於存儲單行中某字符的含量
while(my $seq = <>){
        chomp($seq);
        if($seq =~ m/^>(.*)/){
            $id = $1;
            next;
        }
        @num = ($seq =~ m/(G|C)/g);
        $GC_content{$id} += @num; 
        @num = ($seq =~ m/(.)/g);
        $sequences{$id} += @num; 
}

    foreach(keys(%GC_content)){
        $GC_content{$_} /= $sequences{$_};
    }
my @sort = sort{$GC_content{$b} <=> $GC_content{$a}} keys(%GC_content);
printf("%s\n%.6f%\n", $sort[0], $GC_content{$sort[0]}*100);

 

三、技巧

神奇的perl，神奇的sort！！

對數組（或者哈希）排序獲得下標的方式：

# 數字排序：
my @arr = qw(2 3 41 2 34 ); my @result1 = sort{$a <=> $b} @arr; # 獲得下標：
my @result2 = sort{$arr[$a] <=> $arr[$b]} 0..$#arr; # 獲得key：
my %hash = ( one =>1, two =>5, tree=>9 ); my @result3 = sort{$hash{$a} <=> $hash{$b}} keys(%hash); print "數字排序：@result1\n獲得下標：@result2\n獲得key：@result3\n";

 

 

備注：貼一個感覺不錯的代碼（學習學習）

$/ = '>';
<>; # 讀一次">"前的序列，以免下面代碼出錯
while (<>) {
    chomp;
    my ($id, @ary) = split '\n';
    my $seq = join '', @ary; 
    my $ratio = &GC_content($seq);
    if ($ratio > $highest) {
        $highest = $ratio;
        @result = ($id, $ratio);
    }
}
print join "\n", @result;

sub GC_content {
    my ($seq) = @_;
    my $ratio = $seq =~ s/([CG])/$1/g / length($seq) * 100;
    return $ratio
}