gatk-3.8/R/phasing/scripts/plot.testReadLengths.R

25 lines
783 B
R

theta = 10^-3
Fm_BASE = 392 - 2 * 101 # The mean insert size == 190
Fs = 44
meanDepth = 65
nReadsToPhase = 1
params = paste("meanDepth= ", meanDepth, ", nReadsToPhase= ", nReadsToPhase, ", theta= ", theta, ", Fm_BASE= ", Fm_BASE, ", Fs= ", Fs, sep="")
READ_LENGTHS = seq(from=30, to=1000, by=10)
NUM_READ_LENGTHS = length(READ_LENGTHS)
pPhaseReadLength = as.vector(matrix(data = -1, nrow = 1, ncol = NUM_READ_LENGTHS))
for (i in 1:NUM_READ_LENGTHS) {
Fm = Fm_BASE + 2 * READ_LENGTHS[i]
pPhaseReadLength[i] = pDirectlyPhaseHetPair(meanDepth, nReadsToPhase, READ_LENGTHS[i], theta, Fm, Fs)
}
scatter(READ_LENGTHS, pPhaseReadLength, "testReadLengths", xlab="Read length", ylab="Phaseability", main=params, type="b")
save(list = ls(all=TRUE), file = "testReadLengths.RData")