Aligning short sequences to fastq

Question

[BPors]

Hi,

I am trying to search for the presence of couple sequences (around 400) each with a size of 23 bps,in different fastq files, while allowing 1-2 mismatches at maximum. I am not sure if turning the fastq to a genome(transcriptome) would be a nice approach? I have tried making the fastq -> fasta -> building blast database -> running blastn, however it did not run as my query is not only one sequence.

Example part of my query.file :

ATTTTTCTGAAAAACCCCCTACGA

AACAGGAAGTCAAAAAAAGCCAA

AGGATTTTTTTTTTTCTGGGGACA

The output I am aiming to have is, for each read in my query.file, which of these sequences are having 100% (or having 1-2 mismatches) match in fastq file, and possibly where in the fastq file.

I would appreciate your suggestions! Thank you!

[ATpoint]

You could use bowtie instead of blast. Make a fasta from the fastq, build a bowtie index from it, then align the query. Bowtie has an option that controls how many mismatches are allowed in the seed (-n). As the seed (28bp) is longer than your queries, setting the max seed mismatches to 1 or 2 should be sufficient for your goal.

[BPors]

Thank you for your answer. I would like to try, but I have these reads in just text format, therefore I cannot turn it to fastq. I think in Bowtie I have use reads in fastq format

[ATpoint]

No, several formats are accepted:

-q query input files are FASTQ .fq/.fastq (default) |||| -f query input files are (multi-)FASTA .fa/.mfa |||| -r query input files are raw one-sequence-per-line

[BPors]

Thank you! I have eventually used BBDUK but I will give bowtie a try soon with these options. ( -r).

[ATpoint]

I was not aware of that these is a function in BB. This BB stuff is really a jack-of-all-trades.

[Titus]

Hi,

May be you can try to ta align with bwa aln your 23 bps seq against your fastq files as ref after you transformed it as fasta ?

Best

[BPors]

Thank you for your suggestion. Would this work if my reads are in text format?

Answer 1

[Brian Bushnell]

You can grab the fastq sequences containing these 23-mers with BBDuk like this:

bbduk.sh in=file.fastq outm=matched.fastq ref=23mers.fa k=23 hdist=2

"hdist=2" allows 2 mismatches; you can alteratively set that to 1 or 0. This does not tell you where the match is, but you can do that like this:

bbduk.sh in=matched.fastq out=masked.fastq ref=23mers.fa k=23 hdist=2 kmask=lc

That will convert the matched regions to lowercase.

[BPors]

Thank you! That worked well for me!