Tool:Art: Simulation Tools To Generate Synthetic Next-Generation Sequencing Reads.

Question

[Istvan Albert]

ART is a set of simulation tools to generate synthetic next-generation sequencing reads.

ART simulates sequencing reads by mimicking real sequencing process with empirical error models or quality profiles summarized from large recalibrated sequencing data. ART can also simulate reads using user own read error model or quality profiles. ART supports simulation of single-end, paired-end/mate-pair reads of three major commercial next-generation sequencing platforms: Illumina's Solexa, Roche's 454 and Applied Biosystems' SOLiD. ART can be used to test or benchmark a variety of method or tools for next-generation sequencing data analysis, including read alignment, de novo assembly, SNP and structure variation discovery. ART was used as a primary tool for the simulation study of the 1000 Genomes Project .

ART is implemented in C++ with optimized algorithms and is highly efficient in read simulation. ART outputs reads in the FASTQ format, and alignments in the ALN format. ART can also generate alignments in the SAM alignment or UCSC BED file format.

[Jonathan King]

I'm interested in using ART to simulate Illumina reads at 300bp lengths. This doesn't seem immediately possible given the provided quality profiles. Any idea if it is possible to use ART to simulate 300bp Illumina reads or should I look towards another application?

[Jonathan King]

To answer my earlier question, it is possible to simulate Illumina reads at longer lengths provided you generate the error profiles with ART using input files that have the desired length.

[kshitijtayal]

I am able to synthetic generate reads by ART but i am unable to increase or decrease the error rate in the reads. My question is can we alter the error rate of the synthesis reads using ART. IF yes then how.?

[Al-Hussein]

Istvan Albert May I please have your input/help on an issue that I am facing with ART simulator. I am trying to simulate whole genomes using ART simulator; I am using the HiSeq 2500 model, this is the code I am using:

 ./art_illumina -ss HS25 -sam -i label_Achromobacter.fa -p -l 150 -f 20 -m 400 -s 50 -o paired_data1

However, I keep getting this error:

Error: the number of bases is not equal to the number of quality scores!
qual size: 150,  read len: 143

So, I tried to change the value of the -f quality score, but i still get the same error. May you please help me in figuring out how to overcome/resolve this error? Or point me in a direction that is helpful I am very grateful and appreciative for any input/feedback you would give. Thank you!

[Istvan Albert]

I cannot reproduce this error, the command (other than your reference file) runs fine for me (on a MacOS),

perhaps your input fasta file is not quite right

try it with a different genome to troubleshoot

[Al-Hussein]

Istvan Albert Thank you very much for trying out the command and for your willingness to help. I have going at this for quite sometime now trying to troubleshoot this issue I am facing. However, so far I have not been able to. Please find linkeda fasta file that I am trying. I did like you suggested and tried out different fasta files and see if the error persists. It does. I even tried out some published fasta files and their corresponding code for running the art simulator. However, I also get this error. Which is very mind boggling at the moment.

Linked is one of my fasta files that I am trying to run and this is the code I am using. (The value of the -f I calculated as such: (150*5000000)/3219617 (total genome size) art_illumina --noALN -ss HS25 -i MGYG000000002.fna -p -o art_MGYG00000001. -l 150 -f 169.1815 -m 200 -s 10

I tried it with 10,20,30,40,50 standard deviation value, I tried out various iterations for all values, -s -m -f but I still cant get it to workout. When I run it, on any file it gives me this output Error: the number of bases is not equal to the number of quality scores! qual size: 150, read len: 149

Any help or input is super appreciated and I am very grateful for your time, efforts, and support! Thank you very much!

[Istvan Albert]

I have downloaded the data you mentioned, I have downloaded the latest art binaries for my MacOS and also on a Linux computer, on both when I run it I get:

art/art_illumina -ss HS25 -sam -i MGYG000000002.fna -p -l 150 -f 20 -m 400 -s 50 -o paired_data1

prints:

    ====================ART====================
             ART_Illumina (2008-2016)          
          Q Version 2.5.8 (June 6, 2016)       
     Contact: Weichun Huang <[e-mail]> 
    -------------------------------------------

                  Paired-end sequencing simulation

Total CPU time used: 22.2452

The random seed for the run: 1647370983

Parameters used during run
    Read Length:    150
    Genome masking 'N' cutoff frequency:    1 in 150
    Fold Coverage:            20X
    Mean Fragment Length:     400
    Standard Deviation:       50
    Profile Type:             Combined
    ID Tag:                   

Quality Profile(s)
    First Read:   HiSeq 2500 Length 150 R1 (built-in profile) 
    First Read:   HiSeq 2500 Length 150 R2 (built-in profile) 

Output files

  FASTQ Sequence Files:
     the 1st reads: paired_data11.fq
     the 2nd reads: paired_data12.fq

  ALN Alignment Files:
     the 1st reads: paired_data11.aln
     the 2nd reads: paired_data12.aln

  SAM Alignment File:
    paired_data1.sam

[Istvan Albert]

I also ran:

art/art_illumina --noALN -ss HS25 -i MGYG000000002.fna -p -o art_MGYG00000001. -l 150 -f 169.1815 -m 200 -s 10

and that works as well.

just takes longer

that being said not sure why you set the fold coverage to 169.1815 that does not seem to make much sense.

[Al-Hussein]

Thank you very much for your input and for running the file and testing it out. I am running it on windows. Maybe the problem is due to operating system. I will try to run it on a linux system and see how it goes.

Thank you very much for commenting on the fold coverage, I also think it is a very high number, I am using an equation that uses a number of reads of 5000000. However, I think it is very high and outputs files that are 800mb per file, which is super big.

Can you give a recommendation for fold coverage value, or how to calculate it? I decreased that value down to 1000000 number of reads, however, it is still giving 90-100x coverage. What would you recommend for coverage or a method of calculating the coverage for each genome? (so far, I was under the impression that we need to calculate coverage of each file that we want to simulate using this equation -f specifies the sequencing depth. depth = read length x number of reads / genome size for read length I am using 150 bp x 1000000 number of reads / each individual genome size for each file.

The files I am working on are clustered whole genomes of representative bacterial species.

[Istvan Albert]

first you are setting the fold coverage directly, no other computation is needed

the tool will combine the genome lengths, read length and read numbers to generate the number of reads that produce the preselected coverage, there is no need for you to compute anything else

in general, we set fold coverage as an integer and round number, there is no need for it to be 169.1815

[Al-Hussein]

Thank you very much for the clarification. It makes much more sense now, we just set a whole round integer value for the coverage wanted. For this case, I want the coverage to be 2x, then I just set the -f value to 2. That makes much more sense actually than the 169.1815 that ended up giving a 800mb simulated file. Thank you very much for your help and input! I really appreciate them!

[Al-Hussein]

Istvan Albert Based on your previous recommendations, I was able to run ART on WSL-Ubuntu and have adjusted my fold coverage and the command is working fine. However, I am facing an issue with the output of the ART simulator.

This is an example of the output of one of the genome files I am simulating:

This is an art-simulated genome from a publicly available dataset of MAGs :

As you can see, the file I generated has symbols, letters, and + signs in the middle of the reads. I tried to change the coverage values and standard deviation but I still have these symbols on the output.

this is the code I am using atm: art_illumina --noALN -ss HS25 -i MGYG000000001.fna -p -o art-GCA-900066495. -l 150 -f 2.329469623 -m 200 -s 50 -rs 22

The reason the -f value here is not a solid round integer value, is to ensure that all the different length genomes would end up giving me around the same length of reads/genomes I am using this calculation for each genome -f=(150x50,000)/genome size, this results in files all around 7.9kb in size.

I am simulating them using ART simulator, to make them all the same length, with the same sequencing coverage and have them mimic illumina sequence reads.

Thank you very much in advance for your help, input, and support! Thank you very much for the incredible value you add to this website and to people!

[Istvan Albert]

in general in bioinformatics it is rare that you would need to include an image, all that data is just text that you can copy paste a line from your terminal

it is not clear what you mean, just make sure you understand that ART generates FASTQ files and not FASTA files. The FASTQ format does indeed include various other symbols to represent qualities.

https://en.wikipedia.org/wiki/FASTQ_format

as for the coverage, I think you still copmletely misunderstand what it means, as I said before, you do not need to compute anything

you don't need to multiply this or that. the simulator will do that for you. you tell the simulator what coverage you want and it will generate data to match that.

note how the simulator will generate the number of reads for you, that is the variable parameter that you have no control over.

for other simulators, those that ask you to input the read number N now for those you would need to use a formula to figure out what N gives you a given C coverage based on L and G (read length and genome size). But when the simulator asks you to enter the coverage then you don't need to compute anything else.

[Yancey]

How can I use ART simulator to simulate Illumina reads at the same length, with the different sequencing coverage of each reference sequence? Or can it be done by ART simulator?