perl：计算fasta的GC含量 |
您所在的位置：网站首页 › gc含量怎么算 › perl：计算fasta的GC含量 |
perl：计算fasta的GC含量
|
导读
我的小骆驼已经六岁了，这几天才自己学着走路。 一、目的我的fasta文件是下面那个样子，我要统计每条序列的GC碱基与该序列的总碱基比。顺便比较一下每条序列的GC占比与整个fasta文件的GC占比的大小，结果是下下面那个样子。（其实我是在给circos做配置文件，这个以后再说） 下面那个样子： >k141_641 AAGCAGGTTGGAGAGTACAACCCTACACTTCGTATCGTCATATGCGACACTTTGTATCTG CTCGTCAAAGCTATAGGAAATTTCAATGTTTTTTGCCTTAGATGCCGTTTGG >k141_947 CAGACTATGCCATGAAAGTGCGCAATATAGCCGAAGCAGCCGGTGGAGTGATTCTAGCAG GAACAACTGACGGACTTCTAACCTTCCCAAACACATTTGAGCGTCCGGAGGAAATTAAGT TTTACCGCAATAGTCGCGAGGCAGGAGTCAATACCAGC下下面那个样子： k141_201762 0.436060 big k141_6509 0.466622 big k141_7339 0.403511 small 二、perl脚本 vi contig.gc.pl # 创建文件，并编辑 #!/usr/bin/perl use strict; my %GC_content; # id=>GC 哈希存放id和GC数 my %sequences; # id=>sequence 哈希存放id和碱基数 my ($id, $base_sum, $GC_sum); # id, 碱基数，GC数 my @num; # 中间变量，用于存储单行中某字符数 while(my $seq = ) # 一行行读给$seq { chomp($seq); # 去掉字符串末尾的\n if($seq =~ m/^>(.*)/) # 若该行以“>”开头 { $id = $1; # $1就是第一对小括号中的原符号所对应的匹配内容。 next; # ↓ } @num = ($seq =~ m/(G|C)/g); $GC_content{$id} += @num; # 统计该id序列的GC总数 $GC_sum += @num; # 统计该文件的GC总数 @num = ($seq =~ m/(.)/g); $sequences{$id} += @num; # 统计该id序列的碱基总数 $base_sum += @num; # 统计该文件的碱基总数 } foreach(keys(%GC_content)) # 遍历哈希 { if(($GC_content{$_} / $sequences{$_}) >= ($gc_sum / $base_sum)) # 比较某序列GC含量与总文件GC含量的大小 { printf("%s\t%.6f\tmore\n", $_, $GC_content{$_} / $sequences{$_}); # 如果“>=”，打印该序列的：id，GC含量，more } else { printf("%s\t%.6f\tsmall\n", $_, $GC_content{$_} / $sequences{$_}); # 如果“ output.file # 使用标准输出 ... |
【本文地址】
今日新闻 |
推荐新闻 |