This is a test version of Biostars. For the public version, visit https://www.biostars.org.
pysam duplicated sequence error

Hi,

I did an alignment with STAR, I got the sam file which I can open easily from the terminal:

samtools view -S STAR2Aligned.out.sam

Now I have to work with it in python, so I import the pysam module and try to open the file by:

file = pysam.AlignmentFile("/home/lpp/Desktop/Star_Results/STAR2Aligned.out.sam", "r")

It prints the following error:

[W::sam_hdr_parse] duplicated sequence 'NODE_4_length_21_cov_1.000000'
[W::sam_hdr_parse] duplicated sequence 'NODE_18_length_23_cov_1.000000'

I tried to find the error in different forums but the most similar one to my problem has no answer: http://seqanswers.com/forums/showthread.php?t=58219

Does anybody now where this error is coming from?

Thanks in advance!

sam alignment
• 7,605 views • link
updated by Ram Ram • written by xcalle91 • xcalle91

If you do a 

samtools view -H /home/lpp/Desktop/Star_Results/STAR2Aligned.out.sam | grep NODE_4_length_21_cov_1.000000

then do you get more than one line?

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

It says:

[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
• 0 views • link
updated by Ram Ram • written by xcalle91 • xcalle91

Oh, there should be an -S flag given to samtools as well, mea culpa.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

thankss,

Now It says:

[samopen] SAM header is present: 1457 sequences.
@SQ    SN:NODE_4_length_21_cov_1.000000    LN:41
@SQ    SN:NODE_4_length_21_cov_1.000000    LN:41
• 0 views • link
updated by Ram Ram • written by xcalle91 • xcalle91

I should add that it's likely that you have multiple contigs with the same name in your reference genome. This will simply not work, though STAR won't complain (its output, however, will be broken unless the duplicately-named entries also have duplicate sequence...though even then the MAPQ values and such will be wrong).

• 0 views • link
written by Devon Ryan Devon Ryan

1 answer

Samtools is correct, you have a duplicate header. I suspect that you actually have the same sequence twice, but you might want to double check that. You can use the -t option with samtools and give it a two column file with the deduplicated contig names and length and then you'll be able to get a BAM file. Alternatively, remove the duplicate fasta sequences and remap. I suspect that the former method will be both faster and easier.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks Devon Ryan, that was the problem. I check the contigs and there were two duplicates.

• 0 views • link
updated by Ram Ram • written by xcalle91 • xcalle91

Log in to answer this question.

More posts like this