寫在前面

• 最近在ATAC-Seq中call peaks的時(shí)候發(fā)現(xiàn)需要用到基因組的有效大小（effective genome size）
• 

方法一

自己編寫腳本

#coding=utf-8
import sys
aList=[]
fa_file = sys.argv[1]
with open(fa_file,'r') as f:
    for line in f:
        line = line.strip()
        line = line.upper()
        if not line.startswith(">"):
            baseA = line.count("A")
            baseT = line.count("T")
            baseC = line.count("C")
            baseG = line.count("G")
            aList.extend([baseA, baseT, baseC, baseG])
            # print(aList)
    print("effective_genome_size =", sum(aList))

運(yùn)行腳本，腳本后接基因組名

python genomeSize.py m38

運(yùn)行結(jié)果

effective_genome_size = 2652783500

• 需要統(tǒng)計(jì)總?cè)旧w大小，在循環(huán)中加上baseN = line.count("N")，修改aList.extend([baseA, baseT, baseC, baseG, baseN])

方法二（推薦）

使用統(tǒng)計(jì)軟件gnx

#下載與編譯  
git clone https://github.com/mh11/gnx-tools.git
cd gnx-tools
mkdir bin 
javac -d bin/ src/uk/ac/ebi/gnx/* 
# 沒裝ant，請(qǐng)安裝，鏈接：https://downloads.apache.org/ant/binaries/
# wget https://downloads.apache.org/ant/binaries/apache-ant-1.10.10-bin.tar.gz
# tar -zvxf apache-ant-1.10.10-bin.tar.gz
# ant程序 在/apache-ant-1.10.10/bin里面
ant -f package.xml
#使用方法
java -jar gnx.jar 基因組名

軟件使用

java -jar gnx.jar m38

運(yùn)行結(jié)果

Total number of sequences:          66
Total length of sequences:          2730871774 bp
Shortest sequence length :          1976 bp
Longest sequence length  :          195471971 bp
Total number of Ns in sequences:    78088274
N50:    130694993   (9 sequences)   (1442871972 bp combined)

結(jié)果驗(yàn)證，使用 - 查看有效基因組大小

echo "2730871774-78088274"|bc
2652783500

可見，與上面py腳本運(yùn)行結(jié)果一致