SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bbmap callvariants.sh multisample vcf format, wrong number of fields Greg Bioinformatics 5 04-25-2018 01:00 AM
wrong run parameters. How can I change it after 1 cycle? stasi Illumina/Solexa 4 10-13-2014 01:32 AM
Bowtie default parameters ThePresident Bioinformatics 5 01-23-2014 08:49 AM
Wrong of Tophat's bowtie parameters, a bug chenyao Bioinformatics 3 08-23-2011 06:06 AM
bfast default parameters jy123 Bioinformatics 1 06-04-2011 09:42 PM

Reply
 
Thread Tools
Old 04-01-2019, 01:04 AM   #1
aushev
Member
 
Location: Europe

Join Date: Nov 2009
Posts: 21
Default Wrong mapping with BBmap default parameters

I performed mapping for my RNA-seq data (150 bp reads) with BBMap, just the default parameters (the only special parameter set was low memory):
Code:
bbmap.sh ref=gencode.v29.transcripts.fa path=index_bb_lomem in=sampleA_1.fq.gz out=sampleA.bam overwrite=f usemodulo
When I looked in the resulting sam file, I noticed that for many reads I had CIGAR string looking like this:
Quote:
2=4X6=11I1X7=1X3=1X14=1X2=1X1=1X15=1X3=2X9=1X9=2X21=2X23=1X5=
This how the alignment looked (upper line is the read, lower is suggested hit):
CAGAAGTCAAGGAGTAATTTAGGAGTAATTTTGACTTTCAAGTCTTATTACTTAATAAATACATTTCATAAAGCTACAGCTGCCATAGATAGTGATTCCTCTGATGGATCTGGGCAAAGTCAATTGAAAACCTTCTGGAAAGGAGTCACC
..TGGA......************G.......C...C..............T..G.G...............G...CT.........G.........GT.....................CA.......................T.....

It is clearly incorrect mapping: this read should be just assigned to unmapped. From 1414 reads assigned to my gene of interest (attached), 949 had editing distance of 21 and more (34 reads had distance of 40 and more, up to 357!), i.e. obviously didn't belong to that gene. It means that a lot of the counts I get with the default settings are totally wrong. I am wondering if this is normal result for the mapper and what parameters I need to adjust to avoid this issue - `minratio`, `minind`, or something else?
Attached Files
File Type: gz tigd1.sam.gz (148.7 KB, 0 views)
aushev is offline   Reply With Quote
Old 05-08-2019, 06:53 AM   #2
colindaven
Senior Member
 
Location: Germany

Join Date: Oct 2008
Posts: 415
Default

Bit late, but you could do this

nmtag=f Write NM tags.

and then filter by NM = number of mismatches later.

This might be easier:

editfilter=-1 Ban alignments with more than this many edits.

OR also

idfilter=0 Independant of minid; sets exact minimum identity
allowed for alignments to be printed. Range 0 to 1.
(set to 0.95 ?)

subfilter=-1 Ban alignments with more than this many substitutions.
(set to 5 ?)
colindaven is offline   Reply With Quote
Old 05-08-2019, 07:02 AM   #3
aushev
Member
 
Location: Europe

Join Date: Nov 2009
Posts: 21
Default

Colin, thanks a lot for your input!
(Unfortunately for now I had to switch to STAR because it's hard to get help for BBmap questions while STAR's author is very responsive and helpful)
aushev is offline   Reply With Quote
Reply

Tags
bbmap

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 04:14 PM.


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