Genome 540 Homework Assignment 7

(Winter Quarter 2023)

Due Sunday Feb 26, 11:59pm

• In this assignment, you will create a randomized data set corresponding to the data set from HW6 and (using a new scoring scheme which you compute) run the maximal D-segment algorithm on both the randomized and the real data.
  • First, use the read-start counts you found for the two types of segments in HW6 to determine background frequencies (based on the sums of the counts of the "Non-Elevateed CN segments" and the "Elevated CN segments") and target frequencies (based on the "Elevated CN segments" only) for the 4 possible read start values (0, 1, 2, >= 3). Before computing these frequencies, you should first adjust the count for 0 read starts in background segments (only!) by subtracting 8,422,401. This is the number of N's in the reference chromosome sequence, and should be eliminated because read alignments cannot start at an 'N'.
  • Create a scoring scheme (for each count value 0, 1, 2, 3) based on the background and target frequencies, using LLRs with base 2 logarithms. The reason for using such scores rather than the Poisson-based scores used in HW6 is that the count values in the real data do not fit a Poisson very well, probably because of the fact that read libraries typically involve an amplification step that violates the independence assumption necessary for the Poisson.
  • Write a program that uses the background frequencies above to simulate a sequence of read start counts. The length of this sequence should be the total length of the chromosome used in HW6 minus the number of N's (as given above).

    The following pseudocode gives one approach to creating this sequence:

    N = length of sequence to be simulated
    bkgd[r] = frequency of background sites with r read starts (r = 0, 1, 2, 3).
    for each i = 1...N
        x = random number between 0 and 1 (uniform distribution)
        if x < bkgd[0]
            sim_seq[i] = 0
        else if x < bkgd[0] + bkgd[1]
            sim_seq[i] = 1
        else if x < bkgd[0] + bkgd[1] + bkgd[2] 
            sim_seq[i] = 2
        else 
            sim_seq[i] = 3
    

  • Run your maximal D-segment algorithm on the simulated count sequence with S = -D = 5 and the above scoring scheme. Report a list of pairs, giving for each integer score s = 5, ... 30 the number N_seg(s) of D-segments with score >= s.
  • Run your maximal D-segment algorithm on the 'real data' sequence of read starts used in assignment 6 with the above S and D values, scoring scheme, and list output.
  • Generate a list of ratios from the simulated data output, giving for each s = 5, ... 29, the ratio N_seg(s) / N_seg(s + 1) (rounded to 2 decimal places).

    If there is a 0 in the denominator of your ratio (for cases where there were no D-segments for s), print -1.

  • As discussed in lecture, Karlin-Altschul theory predicts that, for LLR scores using logarithmic base b, the number of D-segments with scores >= s should be proportional to b^-s (b to the power -s; this is the reciprocal of the corresponding LR). Since your scores used logarithmic base 2, if N_seg(s1) is the number of D-segments found with score value >= s1, and N_seg(s2) is the number of D-segments found with score value >= s2, then the ratio N_seg(s1)/N_seg(s2) should be approximately equal to 2^(s2 - s1). Consider the following questions:
    • Does this relationship appear to be true for the simulated data?
    • Is it true for the real data?
    • Would you expect it to be true for the real data?
    • What score threshold is a reasonable one to use for the real data, to ensure a very low false positive rate?
  • Your output should include: the background and target frequencies you found in the first step above (rounded to 4 decimal places) and the scoring scheme you derived from them (rounded to 4 decimal places), as well as the lists of pairs for the simulated and real data from running your program and your answers to the above questions. See the template file for formatting details.
• Use the same template structure for your output on the actual file. Please put everything into ONE file - do not send an archive of files or a tar file. After creating a plain text file (NOT a word processing document file) in this format, compress it (using either Unix compress, or gzip -- if you don't have access to either of these programs, let us know), and send it as an attachment to both Phil at phg@uw.edu and Conor at concamp@uw.edu.