Basic question mapping read/contigs (Bowtie2)

[Alun3.1]

Hi,

I am learning bioinformatics and have a basic question about bowtie2.

I got 2 different sets of single-end sequencing reads:
- set #1: 10,000,000 reads
- set #2: 12,000,000 reads

I concatenated these two datasets together (22,000,000 reads) and assembled them with Trinity (100,000 contigs).

Now, I am trying to know which contigs come from which set of data.

For that I aligned with bowtie2 the indexed contigs to the reads of the first set, and then did the same with the second set separately:
- indexed contigs VS set #1
- indexed contigs VS set #2

I used the bowtie2 --end_to_end option and the --al option in order to output the contig sequences that aligned to the reads, and it returned:
- overall alignment rate of contigs VS set #1 = 80% (8,000,000 sequences)
- aligned contig file set #1: 15,000,000 sequences

- overall alignment rate of contigs VS set #2 = 78% (9,360,000 sequences)
- aligned contig file set #2: 19,000,000 sequences

How can I end up with more contigs than Trinity produced and than the number of reads???

I am clearly doing something wrong.
How could I obtain the sequences of the contigs from set #1 and the sequences of the contigs from set #2?
Is the --al the right option?
Should I start digging the SAM file instead?

Thanks in advance for your help!

[Brian Bushnell]

You should align your reads to the contigs, not the other way around. Generally, though, if you assembled the sequences together, most of the contigs will probably come from both datasets, so I don't think you can answer the question you asked. You can, however, find out which contigs are better represented by which dataset.

Once you do the alignment, you can use pileup (from BBTools) to calculate the coverage of each contig from the two sets:

pileup.sh in=mapped.sam out=stats.txt

You can, alternatively, output the coverage directly from BBMap, using the "covstats" flag. BBMap also has a "kfilter" flag which will prevent a read from being mapped to a location unless there are at least k consecutive matching bases. This is designed to ensure that if you assembled with some kmer length K, only reads with at least K consecutive matching kmers will be mapped to that location - therefore, that specific read actually was used to assemble that specific contig.

[Alun3.1]

Thanks Brian !

You should align your reads to the contigs, not the other way around.

It makes more sense indeed.
However, I also just tried samtools view and got the same number of sequences as bowtie calculated

samtools -c -F 4 <input.sam>

The trick is that I think the sequences in the SAM file are the ones of the reads and not the contigs. The other way around should give me the contig sequences though.

But if I align the reads to the contigs as you mentioned should I use local alignment instead of end-to-end?
Also, how could I discard sequences that align with gaps? Do you think setting --gbar <int> to the longest contig value would work? (--gbar <int> disallows gaps within <int> positions of the beginning or end of the read)

Generally, though, if you assembled the sequences together, most of the contigs will probably come from both datasets

I read several times that the more reads you use the better the assembly, but it is another (complex) topic maybe…