Help with Bowtie2

[aakriti]

Hi all!

I am currently working on a project where I am trying to align assembled scaffolds from the Human Microbiome Project to a gene sequence (agaA from C. striatum). I am having trouble with this since every time I try to make an alignment, I get a result like this:

27528 reads; of these:
27528 (100.00%) were unpaired; of these:
27528 (100.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
0 (0.00%) aligned >1 times
0.00% overall alignment rate

I'm not sure if this is an accurate result. Is there a better way to approach this problem? The reason I'm doing this is because I'm trying to find mutational hot spots in the agaA gene (which I know for sure is found on the skin microbiome) and so I thought that I could perhaps use the data available already in the HMP.

Thanks!

[hui_shi]

So basically, none of your reads were mapped according to what you have shown here.

I think you could try to relax the parameters a bit, so see if any of your reads align to the reference genome.

[aakriti]

Thanks for the response! Where can I set the parameters for the alignment? Is it in the Bowtie2_build file?

[Brian Bushnell]

If the query sequences are longer than the reference sequence, that may cause problems. You might be better off shredding the sequences to something like 200bp before mapping.

[aakriti]

Do you mean dividing the query sequences further? I believe this could be the case since the reference sequence is just one gene.

[Brian Bushnell]

I mean breaking the query sequences into pieces each of length ~200bp, so that they will be shorter then the reference sequence.

[GenoMax]

If the scaffolds are much larger than the gene how about doing the search other way round (gene against the scaffolds)?

[aakriti]

OKay I can write a script to try that. Do you think it would be possible to use the query sequences as the reference and try to align my gene of interest to that? Or would that not work (still trying to understand how the indexing works :|)

[aakriti]

Ah! Genomax, you took the words out of my mouth :P

[GenoMax]

Depending on how many mutations you expect you may want to use a non-NGS aligner.

Brian: Would Pacbio aligner from BBMap work?

[dpryan]

Actually, you can also just use local alignment if the query is longer than the reference. That often works.

[aakriti]

Thank you for all the suggestions! I tried the alignment method suggestions you guys sent me (actually still have to look into BBMap), but am getting the same results (no alignment).

Do you guys have any suggestions on loosening the parameters? Will this affect the integrity of the alignment though?

Also, I was wondering if you guys knew much about the HMP database in general. I'm still trying to get my head around what the scaffolds are referring to and what exactly the assembly fasta files contain. The documentation on this seems a bit sparse/doesn't correspond well with the data on the website. Is the chance high that my gene is actually not even present in the sequences assembled? (Though this seems highly unlikely that it wouldn't be present in any of the scaffolds).

Thanks!

[Brian Bushnell]

Depending on how many mutations you expect you may want to use a non-NGS aligner.

Brian: Would Pacbio aligner from BBMap work?

If the scaffolds diverge sufficiently from the gene sequence, then yes, the BBMap Pac Bio version (mapPacBio8k.sh) or Blast would be better, as they are much more tolerant of low-identity alignments than Bowtie2.

[Mark]

Agreed. Just make a blast db out of your scaffolds and blastn them with your gene