导读
我的小骆驼已经六岁了，这几天才自己学着走路。
一、目的
我的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{$_});
# 如果“<”，打印该序列的：id，GC含量，less
}
}
三、运行脚本
perl contig.gc.pl input.file > output.file
# 使用标准输出
