统计各染色体碱基数
1.切割染色体fasta:
awk '/^>/ { gsub(">","",$1)
FILE=$1 ".fa"
print ">" $1 >> FILE
next}
{ print >> FILE }' genome.fa
2.分别统计:
# 用于fasta格式文件的read数、碱基数、最长的read、最短的read及平均read长度
perl -ne 'BEGIN{$min=1e10;$max=0;}next if ($.%2);chomp;$read_count++;$cur_length=length($_);$total_length+=$cur_length;$min=$min>$cur_length?$cur_length:$min;$max=$max<$cur_length?$cur_length:$max;END{print qq{Totally $read_count reads\nTotally $total_length bases\nMAX length is $max bp\nMIN length is $min bp \nMean length is },$total_length/$read_count,qq{ bp\n}}' input.fa
注意:这个有局限性,不能统计N,下面的脚本可以。
grep:碱基数目和GC含量的统计
grep -v '>' input.fa| perl -ne '{$count_A=$count_A+($_=~tr/A//);$count_T=$count_T+($_=~tr/T//);$count_G=$count_G+($_=~tr/G//);$count_C=$count_C+($_=~tr/C//);$count_N=$count_N+($_=~tr/N//)};END{print qq{total count is },$count_A+$count_T+$count_G+$count_C+$count_N, qq{\nGC%:},($count_G+$count_C)/($count_A+$count_T+$count_G+$count_C+$cont_N),qq{\n} }'
total count is 248956422
GC%:0.417242922380087
用samtools统计:
samtools faidx chr1.fa
chr1 248956422 6 60 61
附加:
# 用于fastq格式文件的read数、碱基数、最长的read、最短的read及平均read长度
perl -ne 'BEGIN{$min=1e10;$max=0;}next if ($.%4);chomp;$read_count++;$cur_length=length($_);$total_length+=$cur_length;$min=$min>$cur_length?$cur_length:$min;$max=$max<$cur_length?$cur_length:$max;END{print qq{Totally $read_count reads\nTotally $total_length bases\nMAX length is $max bp\nMIN length is $min bp \nMean length is },$total_length/$read_count,qq{ bp\n}}' input.fq