simulate_fastq_file

Instructions

Implementation Reference

• src/tools/simulateFastq.ts:74-172 (handler)

  Main handler function for simulate_fastq_file tool. Processes input parameters, generates sequencing reads (single-end or paired-end), calculates coverage, applies mutations and sequencing errors, generates quality scores based on the selected model (illumina/454/ion-torrent/pacbio), and formats output as FASTQ or JSON.

  async handler({ 
    referenceSequence,
    readLength,
    coverage,
    readType = "single-end",
    insertSize = 300,
    insertSizeStd = 50,
    errorRate = 0.01,
    qualityModel = "illumina",
    mutationRate = 0.001,
    seed,
    outputFormat = "fastq"
  }: { 
    referenceSequence: string;
    readLength: number;
    coverage: number;
    readType?: string;
    insertSize?: number;
    insertSizeStd?: number;
    errorRate?: number;
    qualityModel?: string;
    mutationRate?: number;
    seed?: number;
    outputFormat?: string;
  }) {
    const generator = new SequenceGenerator(seed);
    const refLength = referenceSequence.length;
    
    // Calculate number of reads needed
    const totalBasesNeeded = refLength * coverage;
    const effectiveReadLength = readType === "paired-end" ? readLength * 2 : readLength;
    const numReads = Math.ceil(totalBasesNeeded / effectiveReadLength);
    
    const reads: any[] = [];
    
    for (let i = 0; i < numReads; i++) {
      if (readType === "single-end") {
        const read = generateSingleEndRead(
          referenceSequence, 
          readLength, 
          errorRate, 
          mutationRate, 
          qualityModel, 
          generator, 
          i + 1
        );
        reads.push(read);
      } else {
        const readPair = generatePairedEndReads(
          referenceSequence,
          readLength,
          insertSize,
          insertSizeStd,
          errorRate,
          mutationRate,
          qualityModel,
          generator,
          i + 1
        );
        reads.push(readPair);
      }
    }
    
    const stats = {
      totalReads: readType === "paired-end" ? reads.length * 2 : reads.length,
      readLength,
      coverage: Math.round((reads.length * effectiveReadLength / refLength) * 100) / 100,
      readType,
      errorRate,
      mutationRate,
      qualityModel,
      seed: seed || "random",
      citation: "Stephens et al. (2016) 'Simulating Next-Generation Sequencing Datasets from Empirical Mutation and Sequencing Models.' PLOS ONE 11(11): e0167047. DOI: 10.1371/journal.pone.0167047"
    };
    
    let output = '';
    if (outputFormat === 'fastq') {
      if (readType === "single-end") {
        output = reads.map(read => formatFastqRecord(read)).join('\n');
      } else {
        output = reads.map(pair => 
          formatFastqRecord(pair.read1) + '\n' + formatFastqRecord(pair.read2)
        ).join('\n');
      }
    } else {
      output = JSON.stringify(reads, null, 2);
    }
    
    return {
      content: [{
        type: "text",
        text: JSON.stringify({
          statistics: stats,
          reads: outputFormat === 'fastq' ? undefined : reads,
          fastqOutput: outputFormat === 'fastq' ? output : undefined
        }, null, 2)
      }]
    };
  }

• src/tools/simulateFastq.ts:5-73 (schema)

  Tool definition and input schema for simulate_fastq_file. Defines all input parameters including referenceSequence, readLength, coverage, readType, insertSize, errorRate, qualityModel, mutationRate, seed, and outputFormat with their types, constraints, and descriptions.

  definition: {
    name: "simulate_fastq_file",
    description: "Simulate FASTQ sequencing reads with realistic quality scores and error models based on NEAT methodology. Implementation inspired by: Stephens et al. (2016) 'Simulating Next-Generation Sequencing Datasets from Empirical Mutation and Sequencing Models.' PLOS ONE 11(11): e0167047. https://doi.org/10.1371/journal.pone.0167047",
    inputSchema: {
      type: "object",
      properties: {
        referenceSequence: {
          type: "string",
          description: "Reference DNA sequence to generate reads from"
        },
        readLength: {
          type: "number",
          description: "Length of each sequencing read",
          minimum: 50,
          maximum: 300
        },
        coverage: {
          type: "number",
          description: "Target sequencing coverage depth",
          minimum: 1,
          maximum: 1000
        },
        readType: {
          type: "string",
          description: "Type of reads to generate",
          enum: ["single-end", "paired-end"]
        },
        insertSize: {
          type: "number",
          description: "Mean insert size for paired-end reads (ignored for single-end)",
          minimum: 100,
          maximum: 1000
        },
        insertSizeStd: {
          type: "number",
          description: "Standard deviation of insert size for paired-end reads",
          minimum: 10,
          maximum: 200
        },
        errorRate: {
          type: "number",
          description: "Base calling error rate (0-1)",
          minimum: 0,
          maximum: 0.1
        },
        qualityModel: {
          type: "string",
          description: "Quality score model to use",
          enum: ["illumina", "454", "ion-torrent", "pacbio"]
        },
        mutationRate: {
          type: "number",
          description: "Rate of true mutations to introduce (0-1)",
          minimum: 0,
          maximum: 0.05
        },
        seed: {
          type: "number",
          description: "Random seed for reproducible results (optional)"
        },
        outputFormat: {
          type: "string",
          description: "Output format for the reads",
          enum: ["fastq", "json"]
        }
      },
      required: ["referenceSequence", "readLength", "coverage"]
    },
  },

• src/server.ts:107-120 (registration)

  Registration of simulate_fastq_file handler in the MCP server's CallToolRequestSchema switch statement. Routes requests to the simulateFastq.handler with type-validated arguments.

  case "simulate_fastq_file":
    return await simulateFastq.handler(args as {
      referenceSequence: string;
      readLength: number;
      coverage: number;
      readType?: string;
      insertSize?: number;
      insertSizeStd?: number;
      errorRate?: number;
      qualityModel?: string;
      mutationRate?: number;
      seed?: number;
      outputFormat?: string;
    });

• src/utils/sequenceUtils.ts:48-176 (helper)

  SequenceGenerator helper class providing random number generation with optional seed. Used by simulate_fastq_file for deterministic and non-deterministic simulations via the rng() method.

  export class SequenceGenerator {
    public rng: () => number;
  
    constructor(seed?: number) {
      if (seed !== undefined) {
        let s = seed;
        this.rng = () => {
          s = Math.sin(s) * 10000;
          return s - Math.floor(s);
        };
      } else {
        this.rng = Math.random;
      }
    }
  
    generateRandomDNA(length: number, gcContent: number = 0.5): string {
      const gcProb = gcContent / 2;
      const atProb = (1 - gcContent) / 2;
      
      let sequence = '';
      for (let i = 0; i < length; i++) {
        const rand = this.rng();
        if (rand < atProb) {
          sequence += 'A';
        } else if (rand < atProb * 2) {
          sequence += 'T';
        } else if (rand < atProb * 2 + gcProb) {
          sequence += 'G';
        } else {
          sequence += 'C';
        }
      }
      return sequence;
    }
  
    generateRandomProtein(length: number): string {
      let sequence = '';
      for (let i = 0; i < length; i++) {
        const index = Math.floor(this.rng() * AMINO_ACIDS.length);
        sequence += AMINO_ACIDS[index];
      }
      return sequence;
    }
  
    mutateDNA(sequence: string, params: MutationParameters): string {
      const {
        substitutionRate = 0.01,
        insertionRate = 0.001,
        deletionRate = 0.001,
        transitionBias = 2.0
      } = params;
  
      let mutated = sequence.split('');
      
      for (let i = 0; i < mutated.length; i++) {
        if (this.rng() < substitutionRate) {
          const currentBase = mutated[i];
          let newBase: string;
          
          if (this.rng() < transitionBias / (transitionBias + 1)) {
            newBase = this.getTransition(currentBase);
          } else {
            newBase = this.getTransversion(currentBase);
          }
          mutated[i] = newBase;
        }
        
        if (this.rng() < insertionRate) {
          const randomBase = DNA_BASES[Math.floor(this.rng() * DNA_BASES.length)];
          mutated.splice(i, 0, randomBase);
          i++;
        }
        
        if (this.rng() < deletionRate && mutated.length > 1) {
          mutated.splice(i, 1);
          i--;
        }
      }
      
      return mutated.join('');
    }
  
    private getTransition(base: string): string {
      switch (base) {
        case 'A': return 'G';
        case 'G': return 'A';
        case 'T': return 'C';
        case 'C': return 'T';
        default: return base;
      }
    }
  
    private getTransversion(base: string): string {
      const transversions: Record<string, string[]> = {
        'A': ['T', 'C'],
        'T': ['A', 'G'],
        'G': ['T', 'C'],
        'C': ['A', 'G']
      };
      const options = transversions[base] || [base];
      return options[Math.floor(this.rng() * options.length)];
    }
  
    evolveSequence(sequence: string, params: EvolutionParameters): string[] {
      let population = Array(params.populationSize).fill(sequence);
      const generations = [sequence];
      
      for (let gen = 0; gen < params.generations; gen++) {
        population = population.map(seq => 
          this.mutateDNA(seq, { substitutionRate: params.mutationRate })
        );
        
        if (params.selectionPressure) {
          population.sort((a, b) => this.calculateFitness(b) - this.calculateFitness(a));
          population = population.slice(0, Math.floor(params.populationSize / 2));
          population = population.concat(population);
        }
        
        generations.push(population[0]);
      }
      
      return generations;
    }
  
    private calculateFitness(sequence: string): number {
      const gcContent = (sequence.match(/[GC]/g) || []).length / sequence.length;
      return Math.abs(gcContent - 0.5);
    }
  }

• src/tools/simulateFastq.ts:175-375 (helper)

  Supporting helper functions for FASTQ simulation: generateSingleEndRead, generatePairedEndReads, applyMutations, applySequencingErrors, generateQualityScore, reverseComplement, and formatFastqRecord. These implement the core simulation logic including read generation, mutation application, error modeling, quality score generation, and output formatting.

  function generateSingleEndRead(
    reference: string,
    readLength: number,
    errorRate: number,
    mutationRate: number,
    qualityModel: string,
    generator: SequenceGenerator,
    readId: number
  ): any {
    // Random start position
    const maxStart = Math.max(0, reference.length - readLength);
    const startPos = Math.floor(generator.rng() * (maxStart + 1));
    const endPos = Math.min(startPos + readLength, reference.length);
    
    let sequence = reference.substring(startPos, endPos);
    
    // Apply mutations first (true biological variants)
    sequence = applyMutations(sequence, mutationRate, generator);
    
    // Apply sequencing errors
    const { errorSequence, qualities } = applySequencingErrors(
      sequence, 
      errorRate, 
      qualityModel, 
      generator
    );
    
    return {
      id: `sim_read_${readId}`,
      sequence: errorSequence,
      qualities,
      startPos,
      endPos: startPos + errorSequence.length,
      strand: "+",
      mutations: sequence !== reference.substring(startPos, endPos),
      errors: errorSequence !== sequence
    };
  }
  
  function generatePairedEndReads(
    reference: string,
    readLength: number,
    insertSize: number,
    insertSizeStd: number,
    errorRate: number,
    mutationRate: number,
    qualityModel: string,
    generator: SequenceGenerator,
    readId: number
  ): any {
    // Generate insert size with normal distribution approximation
    const actualInsertSize = Math.max(
      readLength * 2, 
      Math.round(insertSize + (generator.rng() - 0.5) * 2 * insertSizeStd)
    );
    
    const maxStart = Math.max(0, reference.length - actualInsertSize);
    const fragmentStart = Math.floor(generator.rng() * (maxStart + 1));
    const fragmentEnd = Math.min(fragmentStart + actualInsertSize, reference.length);
    
    // Read 1 (forward)
    const read1Start = fragmentStart;
    const read1End = Math.min(read1Start + readLength, reference.length);
    let read1Seq = reference.substring(read1Start, read1End);
    
    // Read 2 (reverse)
    const read2End = fragmentEnd;
    const read2Start = Math.max(read2End - readLength, 0);
    let read2Seq = reverseComplement(reference.substring(read2Start, read2End));
    
    // Apply mutations and errors to both reads
    read1Seq = applyMutations(read1Seq, mutationRate, generator);
    read2Seq = applyMutations(read2Seq, mutationRate, generator);
    
    const read1Result = applySequencingErrors(read1Seq, errorRate, qualityModel, generator);
    const read2Result = applySequencingErrors(read2Seq, errorRate, qualityModel, generator);
    
    return {
      read1: {
        id: `sim_read_${readId}/1`,
        sequence: read1Result.errorSequence,
        qualities: read1Result.qualities,
        startPos: read1Start,
        endPos: read1End,
        strand: "+",
        insertSize: actualInsertSize
      },
      read2: {
        id: `sim_read_${readId}/2`,
        sequence: read2Result.errorSequence,
        qualities: read2Result.qualities,
        startPos: read2Start,
        endPos: read2End,
        strand: "-",
        insertSize: actualInsertSize
      }
    };
  }
  
  function applyMutations(sequence: string, mutationRate: number, generator: SequenceGenerator): string {
    const bases = ['A', 'T', 'G', 'C'];
    let mutated = sequence.split('');
    
    for (let i = 0; i < mutated.length; i++) {
      if (generator.rng() < mutationRate) {
        const currentBase = mutated[i];
        const otherBases = bases.filter(b => b !== currentBase);
        mutated[i] = otherBases[Math.floor(generator.rng() * otherBases.length)];
      }
    }
    
    return mutated.join('');
  }
  
  function applySequencingErrors(
    sequence: string, 
    errorRate: number, 
    qualityModel: string, 
    generator: SequenceGenerator
  ): { errorSequence: string; qualities: string } {
    const bases = ['A', 'T', 'G', 'C'];
    let errorSequence = sequence.split('');
    let qualities = '';
    
    for (let i = 0; i < sequence.length; i++) {
      let hasError = generator.rng() < errorRate;
      let quality: number;
      
      if (hasError) {
        // Apply error
        const currentBase = errorSequence[i];
        const otherBases = bases.filter(b => b !== currentBase);
        errorSequence[i] = otherBases[Math.floor(generator.rng() * otherBases.length)];
        
        // Lower quality for error positions
        quality = generateQualityScore(qualityModel, generator, true);
      } else {
        // No error, higher quality
        quality = generateQualityScore(qualityModel, generator, false);
      }
      
      qualities += String.fromCharCode(33 + quality);
    }
    
    return {
      errorSequence: errorSequence.join(''),
      qualities
    };
  }
  
  function generateQualityScore(qualityModel: string, generator: SequenceGenerator, hasError: boolean): number {
    let quality: number;
    
    switch (qualityModel) {
      case "illumina":
        if (hasError) {
          quality = Math.floor(generator.rng() * 20) + 2; // Q2-Q22
        } else {
          quality = Math.floor(generator.rng() * 20) + 20; // Q20-Q40
        }
        break;
      case "454":
        if (hasError) {
          quality = Math.floor(generator.rng() * 15) + 2; // Q2-Q17
        } else {
          quality = Math.floor(generator.rng() * 20) + 15; // Q15-Q35
        }
        break;
      case "ion-torrent":
        if (hasError) {
          quality = Math.floor(generator.rng() * 10) + 2; // Q2-Q12
        } else {
          quality = Math.floor(generator.rng() * 15) + 10; // Q10-Q25
        }
        break;
      case "pacbio":
        if (hasError) {
          quality = Math.floor(generator.rng() * 8) + 2; // Q2-Q10
        } else {
          quality = Math.floor(generator.rng() * 10) + 8; // Q8-Q18
        }
        break;
      default:
        quality = hasError ? 10 : 30;
    }
    
    return Math.min(93, Math.max(0, quality)); // Clamp to valid FASTQ range
  }
  
  function reverseComplement(sequence: string): string {
    const complement: Record<string, string> = { 'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G' };
    return sequence
      .split('')
      .reverse()
      .map(base => complement[base] || base)
      .join('');
  }
  
  function formatFastqRecord(read: any): string {
    return `@${read.id}\n${read.sequence}\n+\n${read.qualities}`;
  }