2023-12-19统计全基因组测序深度

import sys

input_file = sys.argv[1]

step_len = 100000

result = {}

chrom_stats = {

    "Chr1":249246671,

    "Chr2":262532350,

    "Chr3":170554821,

    "Chr4":124586965,

    "Chr5":181296669,

    "Chr6":210921191,

    "Chr7":214854432,

    "Chr8":182071095,

    "Chr9":161119024,

    "Chr10":181937317,

    "Chr11":149641259,

    "Chr12":173678962,

    "Chr13":188826820,

    "Chr14":164794864

    }

for chrom in chrom_stats.keys():

    step= 0

    step_start = 0

    result[chrom]={}

    result[chrom][step]=0

    for line in open(input_file):

        if line.startswith(chrom):

            pos = line.split()[1]

            depth = line.split()[2]

            if step_start < int(pos) < (step_start + step_len):

                result[chrom][step] += int(depth)

            if int(pos) >= (step_start + step_len):

                step += 1

                step_start += step_len

                result[chrom][step] = 0

                result[chrom][step] += int(depth)

    for steps,nums in result[chrom].items():

        print(chrom, steps, nums/100000, sep="\t")