Find reads that completely cover short sequences

Question

[valerie]

Hi everyone,

I have WGS with paired-end reads length of 100 and several reference sequences length of ~70. I want to find those reads (not pairs, just single reads) which completely cover short sequences, i.e. those reads that map to short sequences so that whole sequence is inside this read. How can I do that using bowtie2?

Thanks in advance!

[Pierre Lindenbaum]

context is missing : what is the input for "several reference sequences length of ~70"

[valerie]

I have an index of ~100 of short sequences I am mapping to

[Pierre Lindenbaum]

what is an "index of 100 short sequence" ? blast ? fasta ? bed ?...

[valerie]

Sorry, didn't get your question. It's fasta.

[Friederike]

does every sequence get its own "chromosome"? then you could check for reads that only span one "chromosome" (while also checking for the length of the read)

Answer 1

[Jeffin Rockey]

An approach that I would follow is this.

Do the alignment with the aligner of your choice against your short sequences fasta file.

Make a bed file with the Sequence name, 0, sequence length,

bedtools intersect the above bed and the aligned bam with fraction 1 (100 percent coverage of your regions) using -f option.

Wait for some more time for better answers.If nothing comes up, give this a try.

[Friederike]

a nice suggestion because it's quickly implemented. you may want to double-check that you're fine with mismatches etc. also, what is up with the remaining 30 bp of the reads (100 bp reads for 70 bp inserts?)

[valerie]

The main issue is that when I am using bowtie2 in default mode it doesn't map anything at all to sequences that are shorter then reads. When I try local mode, I get something that is difficult to analyse, because there are a lot of very short alignments. I was wondering if there is a simple way.

[Friederike]

I think you need to explain a bit more in detail what exactly it is you're trying to do. I.e., why are your reads longer than the sequences?