一起生信啦啦啦Shell程序编写生物信息学

统计各染色体碱基数

2018-10-17  本文已影响8人  苏牧传媒

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

上一篇下一篇

猜你喜欢

热点阅读