Question: How to get inserted bases using pysam pileup
0
gravatar for va90
5 months ago by
va9010
va9010 wrote:

Hi,
I am using pysam.pileup to get the base per information from a bam file.
I am using the following code:

import pysam
samfile = bamfh=pysam.AlignmentFile(bam_file,'rb')
def get_base_info(samfile,chrom,pos,map_qual):
for pileupcolumn in samfile.pileup( chrom, pos, pos+1, truncate= True):
    pileupcolumn.set_min_base_quality(0)
    if pileupcolumn.pos == pos:
        for pileupread in pileupcolumn.pileups:
            if pileupread.alignment.mapping_quality >= map_qual:
                if pileupread.alignment.is_reverse:
                    strand = "-"
                else:
                    strand = "+"
                read_len= pileupread.alignment.query_alignment_length
                qname= pileupread.alignment.query_name
                tag= pileupread.alignment.get_tag('RG')
                flag= pileupread.alignment.flag

                if pileupread.is_del and not pileupread.is_refskip:

                    print(flag,pileupread.alignment.mapping_quality,qname,"Deletion",pos)
                else:
                    qual= pileupread.alignment.query_qualities[pileupread.query_position]
                    read_base= pileupread.alignment.query_sequence[pileupread.query_position]
                    print(read_len,strand,flag,pileupread.alignment.mapping_quality,read_base,qual,qname,pos,tag,sep='\t')
# for example for cordinate chr1 at position 100000 for reads with mapping quality higher than 20
get_base_info(samfile,chrom,100000,20)

by this code, I can get deletions, matches, and mismatches. Anyone knows how I can get insertions using pileup?

Thanks,
Vahid.

bam pysam alignment next-gen • 244 views
ADD COMMENTlink modified 3 days ago by Lhl730 • written 5 months ago by va9010

Hi @va90, I think you can do

pileupcolumn.get_query_sequences(mark_matches=True, mark_ends=True, add_indels=True)

This will reproduce mpileup format, from which you can infer all types of variants.

Another shortcut is to use pileupread.indel, it returns positive values incidcating insertions and insert lengh after this base and negative values indicating deletion legth after it.

ADD REPLYlink modified 2 days ago • written 3 days ago by Lhl730
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1956 users visited in the last hour