Entering edit mode
7.0 years ago
jjmmii
•
0
As of the latest version of art_illumina, a popular software to simulate reads with realistic Illumina quality scores, it doesn't have a built-in error profile for MiSeq pair end 300bp. So I have generated my own error profile for that, but how do I specify art_illumina to use this error profile? Because the -ss option only recognizes built-in profiles. I tried:
$ art_illumina -ss ./miseq_pe300_error_profile.txt -l 300 -i 99_otus.fasta -f 1 -na -amp -p -o simulated_pe300_gg_99_otus_16S
but it says it doesn't recognize:
ART does not has a built-in profile for the given sequencing system: ./miseq_pe300_error_profile
List of all built-in profiles are: (...)
I thought of hacking the script and add this custom profile as a built-in but there should be a more convenient way to do this?
Have you tried using an older version of the software such as ART-GreatSmokyMountains-04-17-2016 that has a built-in MiSeq profile?
No, in fact the latest version that I'm using does have MiSeq profile but only up to 250bp. I need a MiSeq 300bp profile. I don't think older versions of art_illumina would have the 300bp profile because it's quite new.
Edited to include correct URL - Well, in that case, you might want to look into other tools such as BBMap - https://sourceforge.net/projects/bbmap/
Official BBMap repository is at: https://sourceforge.net/projects/bbmap/
Thanks @genomax! I've updated my previous comment.
Unfortunately bbmap requires either a pre-made seqs.qual file or a constant given quality score. The next step of my pipeline (QIIME 2) cannot process a fastq with a constant quality score. And I don't know how to write a script to generate a pre-made seqs.qual file with realistic error profiles for input pair-end reads that are joined already (so they are of different lengths). I finally resorted to using unprocessed FASTQ files (which means redoing the preprocessing steps that have already been done in the FASTA files). Anyway thanks for attempting to help.
Have you thought about writing to authors of ART to see if they can help you?
Some random thoughts: Is the profile for 300 bp significant different than the one for 250? This profile is going to be subject to type of libraries being sequenced (e.g. diversity, sequencer) so is that going to make a huge difference for whatever you are trying to do. There are plenty of MiSeq datasets now available. Can you use a sampling of them instead of generating simulated data?
Yes, I wrote to them but no reply :(
Sorry I didn't make myself clear in my question. My ultimate goal isn't to simulate reads. I'm trying to "steal" the quals from simulated reads and add them into the FASTA files with real data. You made a good point that I can steal them from real datasets. But I forgot the second problem which is my FASTA files are joined pair-end reads, so they're actually not 300bp but the size of the 16S amplicon. So I need to find a real dataset of joined sequences in FASTQ format and the amplicon size has to match and hope that every size has a match (because joined sequences' lengths vary by a few base pairs) or write a script to add/subtract a few base pairs. This effort becomes greater than just starting from raw FASTQ files...
One party that I didn't email to yet is the provider of these FASTA files if they have the FASTQ files as well :P So I'll write to them now and see how it goes.
I think of the MiSeq 300 profile as an extension of the 250 profile, as in, the first 250 bp have the same profiles, but from 250 to 300 the quality score continues to drop.
Thanks for listening, haha.
There should be plenty of datasets that can be merged in that size range. You may need to look around on EBI-ENA. If you can get the raw fastq files then that is the simplest way.
Ok! Thanks genomax and Sej Modha.