pysam biostar post

README

python fct_depth_biostars.py ex1.bam > depth_fct.txt

satools depth ex1.bam > depth_samtools.txt
 

fct_depth_biostars.py

#!/usr/bin/python
__author__ = "GENTY Fabien"

import os
import sys
import pysam
import time


def get_chromosomes_names(input):

    # opening the bam file with pysam
    bamfile = pysam.AlignmentFile(input, 'rb')
    # query all the names of  the chromosomes in a list
    list_chromosomes = bamfile.references[0:24]
    list_length = bamfile.lengths[0:24]
    bamfile.close()
    return list_chromosomes, list_length


def count_depth(chr_name, size, input):
    """
    Count the depth of the read. For each genomic coordonate return the
    number of reads
    -----
    Parameters :
        chr : (str) name of the chromosome
    -----
    Returns :
        none
    """
    bamfile = pysam.AlignmentFile(input, 'rb')
    for pileupcolumn in bamfile.pileup(chr_name):
        pos = pileupcolumn.reference_pos + 1
        depth = pileupcolumn.nsegments
        print("{}  {}  {}".format('chr1', pos, depth))


if __name__ == "__main__":
    file = sys.argv[1]
    # query all the names of  the chromosomes in a list
    list_chromosomes, list_length = get_chromosomes_names(file)
    # for chr, size in zip(list_chromosomes, list_length):
    count_depth(list_chromosomes[0], list_length[0], file)