0
$\begingroup$

I am new to python and trying to learn. The below is an attempt to filter out secondary reads in a bam using pysam. I ran mark duplicates and got an error on several bam files due to secondary alignments. I added comments to each line to illustrate my thought process. Thank you :).

The specific mark duplicates error:

Value was put into PairInfoMap more than once

Here is the script:

#! /usr/bin/python  ## call python script
import sys ## import python system functions
import pysam  ## import module

bam = pysam.AlignmentFile(*.bam, "rb")  ## open bam and read
click.echo("Reading BAM file")  ## output message
for read in bam:  ## start loop and iterate over each bam
    if read.is_secondary=true  ## not the primary alignment
    if read.has_tag('XS')  ## look for Secondary Alignment Score tag
        XS = read.get_tag('XS')  ## read and store AS value
    if read.has_tag('AS')  ## look for Alignment Score tag
        AS = read.get_tag('AS')  ## read and store value os XS tag  
    if AS < XS  ## Alignment score less then 
        read.write(read)  ## only print primary alignments
bam.close()  ## end processing
$\endgroup$
  • $\begingroup$ Is there an issue with the script? If so what is it? Are you looking for feedback on the script? $\endgroup$ – Bioathlete Apr 13 at 2:34
3
$\begingroup$

Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the read.has_tag('XS') continue block.

I think that if you filter out is_secondary and is_supplementary you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly.

for read in bam.fetch()  
    if read.is_secondary=true  
        continue
    if read.has_tag('XS') <- This block stores XS and then moves on to the next read
        XS = read.get_tag('XS')  
        continue 
    if read.has_tag('AS') 
        AS = read.get_tag('AS')  
    if AS < XS  <- This block compares the AS of the current read to the XS of the previous read that had one.
        read.write(read)  
bam.close()  
$\endgroup$
  • $\begingroup$ I have not tried to executec it as this is my first python script i am looking for suggestions. Thank you :). $\endgroup$ – justaguy Apr 13 at 5:25
  • 1
    $\begingroup$ I would suggest running the script to determine if it is doing what you want. What I pointed out above is def a flaw in the logic. $\endgroup$ – Bioathlete Apr 13 at 15:21
  • $\begingroup$ I ran the script (modified it in the post) and there are many syntax errors... being new to python I am not sure what is best... what do you suggest? Thanks :). $\endgroup$ – justaguy Apr 15 at 13:12
  • $\begingroup$ I would start will researching how to write to a file. read.write(read) is not correct syntax. You will need to open a file and write to that. Go error by error until you have a working script then determine if the out is correct. $\endgroup$ – Bioathlete Apr 15 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.