Question: How to map longer reads against a local reference using Bowtie
0
gravatar for Apex92
8 weeks ago by
Apex9220
Apex9220 wrote:

I have 40 fasta files with a read length of ~150. A specific part of these 150 bp long reads should contain ~21 bp specific sequences. I have created a local reference.fa file by these 21 bp sequences and want to perform alignment using bowtie1/bowtie2 to see how many reads do have these specific 21 bp sequences with 0 mismatches.

I tried using bowtie 1 by these options but the output was zero alignment: bowtie -v 0 ref -a --best --strata -f *.fasta > bowtie_output .

Based on some forums it seems that bowtie can not handle longer reads against a reference with short reads. Could somebody help me to solve this problem?

I would like to go with zero mismatches (-v 0) and -a --best --strata options. If I have to use bowtie2 which options of it can translate the -v 0 -a --best --strata options in bowtie1?

Thank you.

ADD COMMENTlink modified 8 weeks ago • written 8 weeks ago by Apex9220
1

want to perform alignment using bowtie1/bowtie2 to see how many reads do have these specific 21 bp sequences with 0 mismatches.

You should use bbduk.sh from BBMap suite in filtermode for this purpose, if that is all you want to do.

bbduk.sh in=your.fq literal=21_bp_sequence_1,21_bp_sequence2 etc hdist=0 outm=sequences_matching.fq
ADD REPLYlink modified 8 weeks ago • written 8 weeks ago by GenoMax94k

Thank you for your help. Is there a way to perform this by using bowtie1/bowtie2? I have 64 reference reads with the length of 21 bp and 40 fasta files - In the case of bbduk.sh how can I adjust your suggested command-line and get output for all of my fasta files separately? And the last question of me is that is there an option to increase the number of mismatches in bbduk.sh?

ADD REPLYlink written 8 weeks ago by Apex9220
1

Are you interested in knowing how many reads contain each reference you have? If so you can run each sequence pattern on the command line replacing the literal= sequence each time. If you don't actually need to separate the reads then you can eliminate outm= directive and you will get just the statistics.

If you don't care about individual reference matches, you can create a multi-fasta file with all 64 reads in it and then use the following command line to get all sequences that match one or more of these references).

bbduk.sh in=your.fq ref=multi_fasta_ref.fa hdist=0 outm=sequences_matching.fq
ADD REPLYlink modified 8 weeks ago • written 8 weeks ago by GenoMax94k

Yes, that is exactly what. I want to know - [how many reads contain each reference I have]. My reference is a fasta file and I ran this command for i in `ls -1 *.fasta`; do /Users/bbmap/bbduk.sh in=$i ref=edt_AB_seq.fa hdist=0 outm=$i\sequences_matching.fa; done but I have no results. But for example, when I cat all my reads (40 fasta files) and I grep for one of the reference reads I do e.g. 3 hits - but with bbduk.sh all of the output files are empty.

ADD REPLYlink modified 8 weeks ago • written 8 weeks ago by Apex9220

Here are head of reads: ```

MN00153:75:000H37WNG:1:11102:26354:1789 GCATACACAGTAGCAATTTAGTTTTTCTTCAGTTGCTTTATGACTCTGAAACCTAATTAAATGTAACATATTTTAACCAAAATTCATAGATCGGAAGAGCACACGTCTGAACTCCAGTCACATTACTCGATCTCGTATGCCGTCTTCTGCT MN00153:75:000H37WNG:1:11102:14072:2278 TGACTTCTGAAGGTAATTGGTTGCATCAGATCTTTGTTTTTCAGACCCAGTGGAGATCGGAAGAGCACACGTCTGAACTCCAGTCACATTACTCGATATCGTATGCCGTCTTCTGCTTGAAAAAAAAAATAGGGGGGGGGGGGGGGGGGGG ```

and here is the ref head:

>PLAU_fw
CCAGGTAGTAGTACGTCTGTT
>PU_fw_rc
AACAGACGTACTACTACCTGG
>PU_rv
CGTGCTGCGAGAGTATTATCT
>PU_rv_rc
AGATAATACTCTCGCAGCACG
>AX_fw
CGGCCTCCTCCAATTAAAGAA
ADD REPLYlink modified 8 weeks ago • written 8 weeks ago by Apex9220

I am not sure what you are doing above in the loop.

Since you need to know

how many reads contain each reference I have

you need to take one pattern e.g. CCAGGTAGTAGTACGTCTGTT and then run it against your fasta files to get the reads that match this pattern.

bbduk.sh in=your.fa literal=CCAGGTAGTAGTACGTCTGTT k=10 outm=PLAU_fw.fa hdist=0
ADD REPLYlink modified 8 weeks ago • written 8 weeks ago by GenoMax94k

Can you explain exactly why aligning the short sequences to the long ones won't work for you?

ADD REPLYlink written 8 weeks ago by swbarnes29.4k
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: 2523 users visited in the last hour
_