2
$\begingroup$

I have a quick question about paired-end short reads. I have multiple genomes that were sequenced with paired-end Illumina NextSeq 200 technology, resulting in two fastq files per sample: 001_R1.fastq.gz and 001_R2.fastq.gz, typical for short read data.

I would like to use a k-mer counting and clustering approach to put these genomes together on a phylogenetic tree. I was able to use both reads in one go to count k-mers but, as you can probably guess, this resulted in a very high frequency of k-mers identified and accounted for.

For this specific use case (not interested in changing my approach), would it make more sense to just use a single file for kmer counting (i.e. only the 001_R1.fastq.gz file)?

EDIT: For more context, I'm trying to do an alignment-free approach to produce a phylogenetic distance comparison of these strains. I was using long reads before and it worked fine but Nanopore reads tend to be very messy. I have short-read data for all of the strains and would like to use that instead but that comes at the price of an overload of data.

I suppose my main concern/question is how do I pre-process these short reads without assembling or aligning to a reference to count kmers? Or is that not possible without at least alignment?

$\endgroup$
5
  • $\begingroup$ Do you mean NextSeq 2000? $\endgroup$
    – gringer
    Sep 13, 2023 at 1:45
  • 1
    $\begingroup$ I would be concerned counting them twice. What you could do is map the paired end reads then count kmers. Or you could stitch the overlapping R1 and R2. Or just do not change what you are doing. If accurate mapping is important, I would suggest to map the reads base don R1 and R2 then count kmers without duplicating the count. $\endgroup$
    – Supertech
    Sep 13, 2023 at 2:05
  • 1
    $\begingroup$ Each of the regions sequenced will be covered by tens or hundreds or even thousands of reads depending on coverage. How are you taking that into account? If you are treating the reads as independent sequences, you will be over counting by orders of magnitude whether you use one file or both. $\endgroup$
    – terdon
    Sep 13, 2023 at 9:30
  • $\begingroup$ Would an alignment be the only way to preprocess these short reads for analysis? Ideally, it would be a completely alignment-free approach but do you mean align to a reference or align the R1.fastq to the R2.fastq? @Supertech $\endgroup$
    – rimo
    Sep 13, 2023 at 15:50
  • $\begingroup$ That's what I'm asking here @terdon ... I honestly don't know how to process short reads for this kind of approach. I know there are a lot of reads here because it's very high coverage but I also know I have been asked to not use an assembly because we want to conserve haplotypes and to not use an alignment to a reference to not bias the seqeunce... $\endgroup$
    – rimo
    Sep 13, 2023 at 15:56

2 Answers 2

3
$\begingroup$

I suspect you won't be able to avoid assembling or aligning. If you are working with strains that don't have their own, well studied reference genome and so want to avoid aligning so as not to bias your data towards whichever strain has been chosen as the reference for your species (please correct me if this isn't what's going on), then you will need to assemble a denovo genome from these reads. At the very least you will need to build the longest consensus contigs you can. Unfortunately, this isn't something I have experience with so I can't help you do this beyond telling you that, as I'm sure you already know, there are tools out there that will try to assemble your reads into contigs.

I will instead try to explain why I think you can't avoid assembly. Short read fastq files, paired or not, contain short sequences of ~150-200 base pairs each. Now, depending on the depth of coverage, each region of your genome will be covered by multiple such short reads. Therefore, the data you have in your fastq files will be extremely repetitive and redundant since each of your genome's regions will be represented by multiple reads, and any attempt to find k-mers will be dogged by this redundancy. This is why you got so many more results than you expected and I don't see how using one file only would help at all.

If you had sequenced using an assay that includes unique molecular identifiers (UMIs), then you would be able to extract non-redundant, representative reads for all regions but I cannot think of a way to do this without UMIs off the top of my head. Not saying it isn't possible, only that you need someone more clever than I.

That said, the best I can come up with is that if your samples were sequenced in the same sequencing run and so can be assumed to have the same sequencing biases, you might be able to build your trees despite the redundancy in the data. My thinking here is that what you are really interested in isn't the k-mers themselves or their level of over-representation in each species, but only the difference in the k-mer content across your samples. If so, given that you sequenced all samples in the same run, they should all have the same level of over-representation and any differences might be enough for you to be able to cluster/build your trees.

This seems worth a try to me, but I haven't worked on something like this in many years so if you go down this route, make sure to double check and ensure I'm not leading you astray!

$\endgroup$
1
$\begingroup$

The problem with using single-end short reads is that they're often not long enough to uniquely match to a database (especially for highly-conserved genes, like the ribosomal RNAs). Whether or not this is a problem for your particular situation is unclear; you haven't provided enough context / information for me to work that out. For example, a single end of a 2x50bp run on a NextSeq 2000 might cause problems (especially for medium to large eukaryote genomes), but a 2x300bp run on a NextSeq 2000 is probably okay.

FWIW, Kraken2 has a paired mode, which combines the results from both read ends from the same template source:

Usage of --paired also affects the --classified-out and --unclassified-out options; users should provide a # character in the filenames provided to those options, which will be replaced by kraken2 with "_1" and "_2" with mates spread across the two files appropriately. For example:

kraken2 --paired --classified-out cseqs#.fq seqs_1.fq seqs_2.fq

See the Wiki / manual here:

https://github.com/DerrickWood/kraken2/wiki/Manual#classification

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.