","",$1) FILE=$1 "....">

统计各染色体碱基数

作者: 苏牧传媒 | 来源:发表于2018-10-17 15:13 被阅读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

    相关文章

      网友评论

        本文标题:统计各染色体碱基数

        本文链接:https://www.haomeiwen.com/subject/yjkrzftx.html