assemblyStats.py 文件源码

def calculateN_50(input_file,ratio):

    lengths=[]
    sums=0

    for seq_record in SeqIO.parse(input_file, "fasta"):
        lengths.append(len(seq_record.seq))

    lengths=sorted(lengths, reverse=True)

    contigsMore1000=0
    totalLengthMore1000=0

    for elem in lengths:
        if elem >=10000:
            contigsMore1000+=1
            totalLengthMore1000+=int(elem)
    N_half=sum(lengths)*ratio

    for i in range(0, len(lengths)):
        sums=sums+lengths[i]
        if sums >= N_half:
            return lengths[i],len(lengths),N_half/ratio,contigsMore1000,totalLengthMore1000