一、关于程序:
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";
时间: 2024-12-15 07:16:24