hg19基因组碱基计数

这篇文章的生物学意义可能没有，主要是为了练习下perl语言的使用，对下载得到的hg19人类基因组中A,T,C,G碱基分别进行计数统计。实际上我更想直接看基因组的GC含量，不过如果能得到各个碱基的数目，那也可以了。
下载地址：
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
解压后，再将染色体merge到一个hg19.fa文件，然后开始进行vim goodchromGC.pl进行编写：

#/usr/bin/perl -w

#count Nucleic acid in chrosome

@chr;

@sequence;

%chrom;

while(<>){

    chomp;
    
    if(/^>/){
        if($chrom{$chr}){

            @leng = countBases($chrom{$chr});
            print "$chr @leng\n";

                }
        
        $chr =$_;
        
        $sequence ="";
        
        }
    
    else{

        $chrom{$chr}.=$_;
}

    }


#print %chrom;


#print keys %chrom;

@leng = countBases($chrom{$chr});

print "$chr @leng\n";





sub countBases{

my($dna)=@_;

my $a=0;

my $t=0;

my $c=0;

my $g=0;

my $e=0;

while($dna=~ /a/ig){$a++}

while($dna=~ /t/ig){$t++}

while($dna=~ /c/ig){$c++}

while($dna=~ /g/ig){$g++}

while($dna=~ /[^actg]/ig){$e++}


my(@result)="$a $t $c $g $e";

return @result;

}

计算hg19.fa 中碱基数目

 time cat hg19.fa|sed -n '/>chrM/,/\<>chrUn_g1000211\>/p'|perl ../temp/goodchromGC.pl >partM.txt

因为我用的腾讯云服务器（1核1G），所以一下子运行不起来全部基因组，就采用sed，分段对hg19.fa基因组进行统计，然后进行merge，比较麻烦。如果服务器性能好，那可以直接这样

time cat hg19.fa|perl ../temp/goodchromGC.pl >result.txt

得到的结果大概是这个样子（依次是A,T,C,G,N）

      1 >chrX 45648952 45772424 29813353 29865831 4170000
      2 >chrY 7667625 7733482 5099171 5153288 33720000
      3 >chr10 38330752 38376915 27308648 27298423 4220009
      4 >chr11 38307244 38317436 27236798 27268038 3877000
      5 >chr11_gl000202_random 9226 8978 11254 10645 0
      6 >chr12 0 0 0 0 0
      7 >chr12 38604831 38624517 26634995 26617050 3370502
      8 >chr13 29336945 29425459 18412698 18414776 19580000
      9 >chr14 25992966 26197495 18027132 18071947 19060000
     10 >chr15 0 0 0 0 0
     11 >chr15 23620876 23597921 17247582 17228387 20836626
     12 >chr16 21724083 21828642 17630040 17701988 11470000
     13 >chr17_ctg5_hap1 429214 432922 355909 362783 100000
     14 >chr17 21159933 21206981 17727956 17700340 3400000
     15 >chr17_gl000203_random 12564 12452 6074 6408 0
     16 >chr17_gl000204_random 19702 17322 22701 21585 0

将结果导出到本地，就可以进行下一步分析

人类染色体中各碱基的占比.png

各个染色体ATCG以及N碱基的占比

最后编辑于：2019.04.25 16:53:35