Question: BWA and SOAP2 cannot map a particular sequence which exists in the genome
1
gravatar for francisckantor
6.3 years ago by
European Union
francisckantor10 wrote:

Hi

I have tried to map this short sequence

>sample_1_x12179
TCCTGTACTGAGCTGCCCCGA

which exists on chr8 twice, using various aligners.

BLAT, BLAST, and Bowtie2 found it. BWA MEM and SOAP2 did not. How come? BWA MEM finds it if I include some extra flanking bases. Is it because the sequence is too short?

BWA MEM command:

bwa mem -t 8 $bwaIndex tmp.fasta

 

SOAP command:

./soap2.20 \
 -a tmp.fasta \
 -o ${goutputfile}.tmp \
 -D /mnt/NGS01/ReferenceDB/mirTools/db/genome/hsa.fa.index \
 -M 0 -r 2 -v 5 -p 8 -l 15

I also tried -M 4 (find best hits), but still no results. If I remove the last two bases from the above sequence, then SOAP2 finds a hit on chr11 that requires 2 mismatches.

 

Thank you.

soap2 bwa alignment • 2.2k views
ADD COMMENTlink modified 6.3 years ago • written 6.3 years ago by francisckantor10
2

You're not using the bwa mem command in your example. And it wouldn't work anyway because bwa mem expects a 70bp minimum query sequence length for aligning. Does running bwa aln without the mem like you've shown work okay?

ADD REPLYlink modified 10 months ago by RamRS28k • written 6.3 years ago by Dan D7.1k

Sorry, I did use 'bwa mem', just accidentally posted the wrong command. I'll try 'bwa aln' too. Thank you for your answer.

ADD REPLYlink written 6.3 years ago by francisckantor10

Can you post your SOAP command as well?

ADD REPLYlink written 6.3 years ago by Dan D7.1k

I've edited my post and added the SOAP2 command.

ADD REPLYlink written 6.3 years ago by francisckantor10
3
gravatar for Dan D
6.3 years ago by
Dan D7.1k
Tennessee
Dan D7.1k wrote:

BWA works for me (but not bwa mem, because it requires 70bp minimum query sequence length):

bwa aln $GRCh37 test.fa > test.sai
bwa samse $GRCh37 test.sai test.fa > test.sam

sample_1_x12179      16      8       41518003        0       21M     *       0 0TCGGGGCAGCTCAGTACAGGA   *       XT:A:R  NM:i:0  X0:i:2  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:21 XA:Z:8,+41517962,21M,0;
_____________________________^_______________________________^^^

The documentation for SOAP2 sucks, but I'm guessing it's an issue with your query read length.

ADD COMMENTlink modified 7 months ago by RamRS28k • written 6.3 years ago by Dan D7.1k
2
gravatar for lh3
6.3 years ago by
lh332k
United States
lh332k wrote:
bwa mem -aT20 $GRC37 test.fa

should give you the result and give you all the perfect matches. However, once there are errors in your short sequences, bwa-aln will have a bigger chance to find the hit.

ADD COMMENTlink modified 7 months ago by RamRS28k • written 6.3 years ago by lh332k

Thank you. It seems to work. I'll read more up on the various parameter settings.

ADD REPLYlink written 6.3 years ago by francisckantor10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1485 users visited in the last hour