SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
extract interchoromosomal pairs from a BAM Tanharsh Bioinformatics 3 04-20-2015 04:13 AM
How to extract read count from bam file lylyf1987 RNA Sequencing 0 04-07-2015 11:19 AM
Extract subalignment from SAM/BAM Mesmer Bioinformatics 4 01-17-2015 03:00 AM
extract FASTA from BAM jwhite Bioinformatics 1 05-15-2014 02:33 PM
Extract base from bam file empyrean Bioinformatics 4 07-03-2012 03:47 AM

Reply
 
Thread Tools
Old 06-04-2018, 05:01 AM   #1
Yexzm
Junior Member
 
Location: Paris

Join Date: Dec 2017
Posts: 5
Default Extract the XS field from bam

Hello everyone,

BWA mem generates for each read an "XS" field (the suboptimal alignment score). When I use samtools view, it's presented this way :
- NS500801:90:HY7JVBGXY:2:21205:8003:11253 147 chrM 958 60 76M = 920 -114 CCCCCTCCCCAATAAAGCTAAAACTCACCTGAGTTGTAAAAAACTCCAGTTGACACAAAATAGACTACGAAAGTGG >>;@[email protected]??=??AAC=???>[email protected]>[email protected][email protected]<>>>>>[email protected]@@?A=B=B>>>><@A=B<=;A>[email protected]=;>= BD:Z:[email protected]@GGIHLLLKKLJCKOJLJJ PG:Z:MarkDuplicates RG:Z:id BI:Z:LLLPTSOOSROPQHOTSQOGGNPPQNQPROLPNMMNNFFFFLNNONPOMLNOMLMNEEKLNMOOPOONMHOROPNN NM:i:0 AS:i:76 XS:i:55

Does anyone know an easy way to extract it ? With R ? I mean I know I could use samtools view + awk but it'll take a long time.

Thanks in advance!
Yexzm is offline   Reply With Quote
Old 06-04-2018, 05:51 AM   #2
lindenb
Senior Member
 
Location: France

Join Date: Apr 2010
Posts: 143
Default

using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

Code:
java -jar dist/bioalcidaejdk.jar -e 'stream().forEach(R->println(R.getAttribute("XS")));'  in.bam
lindenb is offline   Reply With Quote
Old 06-07-2018, 04:12 AM   #3
Yexzm
Junior Member
 
Location: Paris

Join Date: Dec 2017
Posts: 5
Default

Hi lindenb, thank you for your answer,

How can I get the read name too ? I would like to have a table the in the first column the read name, and in the second the XS.
Yexzm is offline   Reply With Quote
Old 06-07-2018, 04:19 AM   #4
lindenb
Senior Member
 
Location: France

Join Date: Apr 2010
Posts: 143
Default

> How can I get the read name too ?

Code:
... printl(R.getReadName()+" "+R.getAttribute("XS")
lindenb is offline   Reply With Quote
Old 06-25-2018, 04:39 AM   #5
Yexzm
Junior Member
 
Location: Paris

Join Date: Dec 2017
Posts: 5
Default

Thank you very much for your help!

The output file is too big, I'm trying to get the chromosome too so that I can separate it per chromosome. I tried "getReferenceIndex" but it returns "null". Do you know how I could do ?
Yexzm is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 11:29 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO