sam file without XS attribute

Question

[blooming.daisy333]

im newbie to NGS and linux. Im using sam file output from MapSplice in Cufflinks and facing the following problem :

**> Processing Locus chr9:43969287-43974380      [***********************  ]  95%SAM error on line 88768364: found spliced alignment without XS attribute
SAM error on line 88768374: found spliced alignment without XS attribute
SAM error on line 88768381: found spliced alignment without XS attribute
SAM error on line 88768383: found spliced alignment without XS attribute
SAM error on line 88768384: found spliced alignment without XS attribute
SAM error on line 88768389: found spliced alignment without XS attribute**

while the sam output file looks like this where some reads have XS attribute while others dont have:

 TTGGCTAATCGTCAACCTGGAAGAGCTTTCCTCT      @CCFDFFFGHHHFIIIIIJIJIIJGGIIJJCGHIJJJHIJJJC)*9BGIJIJIIIGIIIGFGIJIJJIJJJJHGHHFFDFECDDBCCBDDDDCDDDDD      NM:i:0  IH:i:1  HI:i:1  XS:A:+  XF:Z:GTAG,
    SRR959590.257   99      chr8    78240758        113     80M9S   =       78240861        189     TACTGATGATGCCACTTGCCTGGCATGAAGAGGAGATGGTGGAGGTGGAGGAGGTGCAGATTTAGATGCTGAATTTACCTCAAAAGTGG       CCCFFFFFGHHHHJJJJJJJJJJJJJJJJJGIJEHDHII?FHDHI0BG?FHGHI;FHIJJJJJIJHHHHHHFFFFDFFDE(.;CD>-;@       NM:i:1  IH:i:1  HI:i:1
    SRR959590.257   147     chr8    78240861        118     86M     =       78240758        -189    GGCGGGGGTGGTGGGGTTACAGGAGATCCAGCTCTAATAGGTGCAGGAGAAGGGGGAGGTGGAGGTATAGGTGTGGCAGGAGTTGT  0''DDDDDDDCBFHHGHHHEHHGHEFCHDCGC=)FCGDAF?0?0>IHFGEBJJJJIGJJJIJJJJIJJJJJJJHHHHHFFFFFCCC  NM:i:0  IH:i:1  HI:i:1
    SRR959590.377   83      chr13   49046592        139     82M110N19M      =       49046362        -441    ATTTAGCATTTACTGTCAAATATCATAGGGATCCGCAAACTGACTCTGCTAGAATTGTGGGTTTTGAGGTTAAACCATACAGTATTAAACATGAATATGAG   CDDDDDEFDCDCDCDEFFFFDEEFFDFFHHHJIJJIHGDJIHFJJJJJJIJIJJJIJIJJIIJJJJJIJJJGJHJJJGHJJJJJJJJJHHHHHFFFFFCCC   NM:i:0  IH:i:1  HI:i:1  XS:A:+  XF:Z:GTAG,
    SRR959590.377   163     chr13   49046362        135     98M     =       49046592        441     TAACCTTCCCCTTGTTGTCCCCATTAGAAGATTAGATCAGGATTCTTCTACTGTTTATCAGCTTGGCTACCATGTTGGGCTCAAAGGGCAATATAGTG      CCCFFFFFHHHGHIJIJJIJJJIIJGIIGIIIGIJJGIIIJI3B90?BGEBFDGGIGIGEIIIJEIIJJJIIIIIIGFEHFFFD).>7;?BDCDAABC      NM:i:0  IH:i:1  HI:i:1
    SRR959590.491   99      chr10   36684807        121     85M2I14M        =       36684912        155     TGATTCATTATTTTTAGCGTCTCTAGTTCAAAACCGAACATGAAATTTTGGTTTCATTCGGCTCCTTTATGGAAGATGGATAAAATAAAAAAAAAAGAAAT   @@@FDFFFDHGDHICEHIIGDGGHDEHGGGEHGEEG>@D@DFGGGGIIHIIFIGGGHIIID<DA?EHFBBEFC>CACDEDDDDDADCDDDDDDDDD&8??:   NM:i:2  IH:i:4  HI:i:1
    SRR959590.491   355     chr11   34180669        121     85M2I14M        =       34180774        155     TGATTCATTATTTTTAGCGTCTCTAGTTCAAAACCGAACATGAAATTTTGGTTTCATTCGGCTCCTTTATGGAAGATGGATAAAATAAAAAAAAAAGAAAT   @@@FDFFFDHGDHICEHIIGDGGHDEHGGGEHGEEG>@D@DFGGGGIIHIIFIGGGHIIID<DA?EHFBBEFC>CACDEDDDDDADCDDDDDDDDD&8??:   NM:i:2  IH:i:4  HI:i:2
    SRR959590.491   355     chr8    46854583        121     85M2I14M        =       46854688        155     TGATTCATTATTTTTAGCGTCTCTAGTTCAAAACCGAACATGAAATTTTGGTTTCATTCGGCTCCTTTATGGAAGATGGATAAAATAAAAAAAAAAGAAAT   @@@FDFFFDHGDHICEHIIGDGGHDEHGGGEHGEEG>@D@DFGGGGIIHIIFIGGGHIIID<DA?EHFBBEFC>CACDEDDDDDADCDDDDDDDDD&8??:   NM:i:2  IH:i:4  HI:i:3
    SRR959590.491   355     chr8    85697578        121     85M2I14M        =       85697683        155     TGATTCATTATTTTTAGCGTCTCTAGTTCAAAACCGAACATGAAATTTTGGTTTCATTCGGCTCCTTTATGGAAGATGGATAAAATAAAAAAAAAAGAAAT   @@@FDFFFDHGDHICEHIIGDGGHDEHGGGEHGEEG>@D@DFGGGGIIHIIFIGGGHIIID<DA?EHFBBEFC>CACDEDDDDDADCDDDDDDDDD&8??:   NM:i:2  IH:i:4  HI:i:4
    SRR959590.491   147     chr10   36684912        61      48S50M  =       36684807        -155    GGCTCCTTTATGGTAGATGGATAAAATAAAAAAAAAAGAAATAAATAAAATAGCATATAAATAAGATAAAGATAAGGCAGGAACGAAATGAAATAAGA      DDA?5(CCCC@:5(@ADDCCDCCAEC?CFJIJJIHCIEJIIEIHGBBHIGCH??IIIGIIGIJJIJIJJJJGJJJJJJIIGGHJIHHHHHFFFFFCCC      NM:i:1  IH:i:4  HI:i:1
    SRR959590.491   401     chr11   34180774        61      48S50M  =       34180669        -155    GGCTCCTTTATGGTAGATGGATAAAATAAAAAAAAAAGAAATAAATAAAATAGCATATAAATAAGATAAAGATAAGGCAGGAACGAAATGAAATAAGA      DDA?5(CCCC@:5(@ADDCCDCCAEC?CFJIJJIHCIEJIIEIHGBBHIGCH??IIIGIIGIJJIJIJJJJGJJJJJJIIGGHJIHHHHHFFFFFCCC      NM:i:1  IH:i:4  HI:i:2
    SRR959590.491   401     chr8    46854688        61      48S50M  =       46854583        -155    GGCTCCTTTATGGTAGATGGATAAAATAAAAAAAAAAGAAATAAATAAAATAGCATATAAATAAGATAAAGATAAGGCAGGAACGAAATGAAATAAGA      DDA?5(CCCC@:5(@ADDCCDCCAEC?CFJIJJIHCIEJIIEIHGBBHIGCH??IIIGIIGIJJIJIJJJJGJJJJJJIIGGHJIHHHHHFFFFFCCC      NM:i:1  IH:i:4  HI:i:3
    SRR959590.491   401     chr8    85697683        61      48S50M  =       85697578        -155    GGCTCCTTTATGGTAGATGGATAAAATAAAAAAAAAAGAAATAAATAAAATAGCATATAAATAAGATAAAGA

can anyone kindly help me how to sort out this problem?????

[GenoMax]

im newbie to NGS and linux

There is no need to include that line in each thread you create. Even experienced folks need help at times.

BTW: You can subscribe to threads you create by using "by email" option where you see the dropdown for "follow via". Bookmark not needed.

[Kevin Blighe]

...and remember that a 'blooming daisy' will eventually produce beautiful flowers and will stand on its own 2 feet.

[blooming.daisy333]

ooooo thanks for your encouragement. always stay blessed..

[blooming.daisy333]

Thank you so much for the guidance. I ll keep these points in mind.

[blooming.daisy333]

Thank you so much for the guidance. I ll keep these points in mind.

Answer 1

[i.sudbery]

Where your sequencing is unstranded, mappers can still add a strand to a read if it crosses a splice site, and that splice site has a canonical splice site sequence - it can use this sequence to work out which way the RNA is going.

Mappers such as TopHat (which cufflinks was originally designed to work with) only accepted splice sites with one of the canonical splice sites (i'm guessing), thus all spliced reads would contain the XS tag. This is definately true for HiSat when you run it in cufflinks compatible mode.

My guess is that MapSplice is allowing reads without the canonical splice site, and thus there are some reads that it can't orientate, so it isn't adding a tag.

As for how to solve it.... you can either find a way to tell MapSplice not to map non-canonical splice sites, use another mapper or filter out reads with non-canonical splice sites. They are probably not real anyway. There is probably some samtools foo for this, but I would do it in python with pysam:

from pysam import AlignmentFile

inBam = AlignmentFile("mybam.bam")
outBam = AlignmentFile("mybam.filtered.bam", "wb", template=inbam)

for read in inbam:

    # output all unspliced reads
    if 'N' not in read.cigarstring:
        outBam.write(read)

   # read is spliced check for XS tag
   if read.has_tag('XS'):
       outBam.write(read)

outBam.close()

[blooming.daisy333]

thank you so much for the kind help. i wonder when the filter for non canonical sites was off in MapSplice then why it has added those reads ??? I would try your script and will discus the output.