1. Modify this code so that it computes the score for each alignment using the f
ID: 3821898 • Letter: 1
Question
1. Modify this code so that it computes the score for each alignment using the following scoring:
match = +1 mismatch = 0 gap = -1
The program should find the optimal alignment(s) and print them with their score.
# Program to compute and print a dot plot of two nucleotide sequences
# Function to initialize a dot plot of a certain size
def initDotPlot(rows, cols):
# Initialize the dot_plot matrix
empty_dot_plot = []
# For each row in the matrix
for i in range(rows):
# Initialize the row
row = []
# For each col in the matrix
for j in range(cols):
# Insert a blank in the row
row.append("-")
empty_dot_plot.append(row)
return empty_dot_plot
# Function to fill the dot plot
def computeDotPlot(sideStr, topStr, dot_plot):
# For each character in str2 (each row)
for i in range(len(sideStr)):
# For each character in str1 (each col)
for j in range(len(topStr)):
# If the character in str1 matches the character in str2
if topStr[j] == sideStr[i]:
# Change the value in that position to a *
dot_plot[i][j] = ' * '
return dot_plot
# Function to print the dot plot of two sequences
def printDotPlot(sideStr, topStr, dot_plot):
row = ""
for ch in topStr:
row = row + " " + ch + " "
print(" ", row)
row = ""
for i in range(len(dot_plot)):
row = row + sideStr[i] + " "
for j in range(len(dot_plot[i])):
row = row + " " + dot_plot[i][j] + " "
print(row)
row = ""
# Main function
def main():
topStr = input("Enter top string: ")
sideStr = input("Enter side string: ")
empty_dot_plot = initDotPlot(len(sideStr), len(topStr))
dot_plot = computeDotPlot(sideStr, topStr, empty_dot_plot)
printDotPlot(sideStr, topStr, dot_plot)
.... Then, add two counters and some print statements to your code to count how much work is done by this algorithm to find the optimal alignment. The main units of work in this program are adding a gap to a sequence, and comparisons to compute the score. You may use global variables to store the amount of work that is done. Test your result on the following sequences:
AAGGTAGCCTAACGTCCACTTTACCC
AGTAAGGTACCTACCTCAACTTCA
How much work does the algorithm do on these sequences?
NOTE: To use a global variable to represent the amount of work required to insert gaps into the sequences, include the following line of code in each function where you reference the variable:
global gap_work
Similarly for the variable to count the amount of work required to compare the characters in the sequences.
global compare_work
The following is an example of what your output might look like for several test cases:
main Enter string 1 GGAGTATAA Enter string 2: GCGTTA There are 84 alignments The following alignments are optimal GGAGTATAA -GCGT-T-A. Score 2 GGAGTATAA -GCGT-TA- Score 2 GGAGTATAA. G-CGT-T-A Score 2 GGAGTATAA G-CGT-TA- Score 2 GGAGTATAA. GC-GT-T-A Score 2 GGAGTATAA. GC-GT-TA- Score 2 Amount of work done Cgaps 279 Amount of work done (comparisonsExplanation / Answer
setClass(
Class = "C1.dotplot",
contains = "ADEg.C1"
)
setMethod(
f = "initialize",
signature = "C1.dotplot",
definition = function(.Object, data = list(score = NULL, at = NULL, frame = 0, storeData = TRUE), ...) {
.Object <- callNextMethod(.Object, data = data, ...) ## ADEg.C1 initialize
.Object@data$at <- data$at
validObject(.Object)
return(.Object)
})
setMethod(
f = "prepare",
signature = "C1.dotplot",
definition = function(object) {
nameobj <- deparse(substitute(object))
## pre-management of graphics parameters
oldparamadeg <- adegpar()
on.exit(adegpar(oldparamadeg))
adegtot <- adegpar(object@adeg.par)
if(object@data$storeData) {
score <- object@data$score
at <- object@data$at
} else {
score <- eval(object@data$score, envir = sys.frame(object@data$frame))
at <- eval(object@data$at, envir = sys.frame(object@data$frame))
}
score <- as.matrix(score)[, 1] ## to manage 'score' when it is a data.frame with only one column
## change some defaults
adegtot$p1d$rug$draw <- FALSE
## object modification before calling inherited method
object@adeg.par <- adegtot
callNextMethod() ## prepare graph
if(object@adeg.par$p1d$horizontal && is.null(object@g.args$ylim))
object@g.args$ylim <- setlimits1D(min(at), max(at), 0, FALSE)
if(!object@adeg.par$p1d$horizontal && is.null(object@g.args$xlim))
object@g.args$xlim <- setlimits1D(min(at), max(at), 0, FALSE)
assign(nameobj, object, envir = parent.frame())
})
setMethod(
f = "panel",
signature = "C1.dotplot",
definition = function(object, x, y) {
## Drawing dotchart
## x is the index
## y is the score
## get some parameters
pscore <- object@adeg.par$p1d
ppoints <- lapply(object@adeg.par$ppoints, FUN = function(x) {rep(x, length.out = length(x))})
plines <- lapply(object@adeg.par$plines, FUN = function(x) {rep(x, length.out = length(x))})
## reorder the values
y <- y[order(x)]
x <- sort(x)
## Starts the display
## depends on the parametres horizontal
## rug.draw and reverse are always considered as FALSE
if(pscore$horizontal) {
x.tmp <- y
y.tmp <- x
panel.segments(object@adeg.par$porigin$origin[1], y.tmp, x.tmp, y.tmp, lwd = plines$lwd, lty = plines$lty, col = plines$col)
} else {
x.tmp <- x
y.tmp <- y
panel.segments(x.tmp, object@adeg.par$porigin$origin[1], x.tmp, y.tmp, lwd = plines$lwd, lty = plines$lty, col = plines$col)
}
panel.dotplot(x = x.tmp, y = y.tmp, horizontal = pscore$horizontal, pch = ppoints$pch, cex = ppoints$cex, col = ppoints$col, alpha = ppoints$alpha, col.line = "transparent")
})
s1d.dotplot <- function(score, at = 1:NROW(score), facets = NULL, plot = TRUE, storeData = TRUE, add = FALSE, pos = -1, ...) {
## evaluation of some parameters
thecall <- .expand.call(match.call())
score <- eval(thecall$score, envir = sys.frame(sys.nframe() + pos))
## parameters sorted
sortparameters <- sortparamADEg(...)
## facets
if(!is.null(facets)) {
if(NCOL(score) == 1)
object <- multi.facets.C1(thecall, sortparameters$adepar, samelimits = sortparameters$g.args$samelimits)
else
stop("Facets are not allowed with multiple scores")
}
## multiple scores
else if(NCOL(score) > 1) {
object <- multi.score.C1(thecall)
}
## simple ADEg graphic
else {
if(length(sortparameters$rest))
warning(c("Unused parameters: ", paste(unique(names(sortparameters$rest)), " ", sep = "")), call. = FALSE)
## creation of the ADEg object
if(storeData)
tmp_data <- list(score = score, at = at, frame = sys.nframe() + pos, storeData = storeData)
else
tmp_data <- list(score = thecall$score, at = thecall$at, frame = sys.nframe() + pos, storeData = storeData)
object <- new(Class = "C1.dotplot", data = tmp_data, adeg.par = sortparameters$adepar, trellis.par = sortparameters$trellis, g.args = sortparameters$g.args, Call = match.call())
## preparation
prepare(object)
setlatticecall(object)
if(add)
object <- add.ADEg(object)
}
if(!add & plot)
print(object)
invisible(object)
}
Related Questions
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.