simulate_ld_motif.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import StatsBase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this function takes the length as an argument to generate sequences\n",
    "# note, this can also be used to generate a motif\n",
    "function GenerateSequence(Length)\n",
    "    return collect(Random.randstring(\"ACGT\", Length))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create sequences as a vector of char\n",
    "# GenerateSequence(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a vector of vector{chars}, each with the same length\n",
    "function GenerateSequences(NumberOfSequences, Length)\n",
    "    return Sequences = map(i -> GenerateSequence(Length), 1:NumberOfSequences)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test it out\n",
    "GenerateSequences(2, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function Mutate!(Sequence, Distance)\n",
    "    # return the sequence if the Distance is 0\n",
    "    if Distance == 0\n",
    "        return Sequence\n",
    "    end\n",
    "    # choose the sites to mutate\n",
    "    posToMutate = StatsBase.sample(1:length(Sequence), Distance, replace = false)\n",
    "    \n",
    "    # pick the letters that can be used at each position\n",
    "    basesAtPositions = [string(i) for i in Sequence[posToMutate]]\n",
    "\n",
    "    lettersToFill = [\"ACGT\" for i in 1:length(Sequence)]\n",
    "    \n",
    "    # figure out what letters each positions can be changed to\n",
    "    lettersToFill = map((x, y) -> replace(x, y => \"\"), lettersToFill, basesAtPositions)\n",
    "    \n",
    "    # select one character from each of these positions and place them where they should be\n",
    "    # use only to cast the string to a character\n",
    "    # https://stackoverflow.com/questions/59946081/julia-convert-string-to-char-or-convert-arraysubstringstring-1-to-char\n",
    "    Sequence[posToMutate] = map(x -> only(Random.randstring(x, 1)), lettersToFill)\n",
    "    \n",
    "    return Sequence\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate test data\n",
    "# distance is the number of mutations in each implanted motif\n",
    "function GenerateTestData_ld(NumberOfSequences, LengthMotif, LengthSequences, Distance)\n",
    "    # do some error checking to make sure the values provided are valid\n",
    "    # specifically, Disance <= LengthMotif\n",
    "    if Distance > LengthMotif\n",
    "        error(\"The Distance if larger than the LengthMotif.\")\n",
    "    end\n",
    "\n",
    "    # LengthMotif <= LengthSequences\n",
    "    if LengthMotif > LengthSequences\n",
    "        error(\"The LengthMotif is longer than the LengthSequences\")\n",
    "    end\n",
    "    \n",
    "    # NumberOfSequences >= 2\n",
    "    if NumberOfSequences <= 1\n",
    "        error(\"The NumberOfSequences is 1, which is too small for motif detection\")\n",
    "    end\n",
    "\n",
    "    # generate the input sequences\n",
    "    sequences = GenerateSequences(NumberOfSequences, LengthSequences)\n",
    "    \n",
    "    # make the motif\n",
    "    motif = GenerateSequence(LengthMotif)\n",
    "    \n",
    "    # generate a mutated motif for each sequence to implant\n",
    "    motifs = [copy(motif) for i in 1:NumberOfSequences]\n",
    "    \n",
    "    # broadcast the mutated motif\n",
    "    motifs .= Mutate!.(motifs, Distance)\n",
    "    \n",
    "    # now implant the mutated motifs and record where we implant it\n",
    "    motifStarts = rand(1:(LengthSequences - LengthMotif + 1), NumberOfSequences)\n",
    "    \n",
    "    # now place the motifs into the sequences vector\n",
    "    for i in 1:length(sequences)\n",
    "        sequences[i][motifStarts[i]:(motifStarts[i] + LengthMotif - 1)] = motifs[i]\n",
    "    end\n",
    "    \n",
    "    # now return all the things we may want later\n",
    "    return (motif = motif, motifs_starts = motifStarts, motifs_implanted = motifs, sequences = sequences)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GenerateTestData_ld(10, 5, 10, 2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}