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会出现计算错误，知道的朋友欢迎留言）

 1 use strict;

 2 my %GC_content; # id=>GC_content

 3 my %sequences; # id=>sequence

 4 my ($id, $sum); # id, 每个序列的字符个数

 5 my @num; # 中间变量，用于存储单行中某字符的含量

 6 while(my $seq = <>){

 7         chomp($seq);

 8         if($seq =~ m/^>(.*)/){

 9             $id = $1;

10             next;

11         }

12         @num = ($seq =~ m/(G|C)/g);

13         $GC_content{$id} += @num;

14         @num = ($seq =~ m/(.)/g);

15         $sequences{$id} += @num;

16 }

17

18     foreach(keys(%GC_content)){

19         $GC_content{$_} /= $sequences{$_};

20     }

21 my @sort = sort{$GC_content{$b} <=> $GC_content{$a}} keys(%GC_content);

22 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

}

巴特西

perl练习——FASTA格式文件中序列GC含量计算&perl数组排序如何获得下标或者键

最新文章

热门文章