Entering edit mode
5.5 years ago
va90
▴
50
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.
Hi @va90, I think you can do
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.