Sid vs Chia, Round I

HapMap2 10 lanes from 10 separate lines, checking for UniqueCRM2 abundance with Chia et al. (2012) estimates.

Get raw data, add to dataframe.

setwd("/Users/Sid/Projects/")
# l = the inbred maize lines corresponding to the the reads.
l = c("CAUZHENG58", "CAU5003", "CAU178", "CAUMO17", "CAUCHANG72", "CAU478", 
    "B73", "M37W", "P39", "NC358")
# rb = reads mapped per lane with BWA
rb = c(1993, 6850, 606, 5935, 13708, 738, 17819, 20, 9024, 25)
# rm2 = reads mapped per lane with Mosaik2 (These reads were mapped using
# Jose's parameters)
rm2 = c(0, 1712, 140, 301, 587, 35, 19111, 21, 9729, 17)
# rm2d = reads mapped per lane with Mosaik2 default.
rm2d = c(0, 121, 10, 6, 3, 0, 16105, 16, 7825, 13)
# trb = total reads BWA
trb = c(1404397, 7483130, 710098, 4884680, 7832364, 494090, 17608448, 20016, 
    16964382, 26880)
# trm2 = total reads Mosaik 2.
trm2 = c(1403910, 7483130, 710098, 4882699, 7828627, 493928, 17608417, 20016, 
    16964376, 26880)
# t = trancript (e.g Unique CRM2, 3068 bp)
t = 3068
# tpkb = trancript length / kilo base
tpkb = (t/1000)
# trbpM = total reads BWA / million
trbpM = c(trb/1e+06)
# trmpM = total reads Mosaik2 / million
trmpM = c(trm2/1e+06)
# chia = Chia RPKM values.
chia = c(1069.9625, 1008.2056, 1056.06, 910.1055, 1276.0608, 1099.9699, 1112.4219, 
    1077.9139, 661.5023, 1115.9119)
# df is a dataframe of mapped reads from each aligner.
df = data.frame(l, rb, rm2, rm2d, trb, trm2, chia)

Calculate RPKM and compare w/ chia

RPKM = reads per kilo base per million

R² values for each plot.

bwa_rpkm = df$rb/tpkb/trbpM
cor(bwa_rpkm, chia)

## [1] 0.7031

mosaik_rpkm = df$rm2/tpkb/trmpM
cor(mosaik_rpkm, chia)

## [1] -0.06028

cor(bwa_rpkm, mosaik_rpkm)

## [1] -0.5177

mosaik_def_rpkm = df$rm2d/tpkb/trmpM
cor(mosaik_def_rpkm, chia)

## [1] -0.07266

cor(mosaik_def_rpkm, mosaik_rpkm)

## [1] 0.9875

cor(bwa_rpkm, mosaik_def_rpkm)

## [1] -0.4619

Plot: Chia vs BWA

plot(chia ~ bwa_rpkm, main = "Chia vs BWA", xlab = "Sid BWA", ylab = "Chia", 
    pch = 16)
legend("bottomright", "R^2=0.6612", cex = 0.75)
abline(lm(chia ~ bwa_rpkm))

Plot: Chia vs Mosaik2

plot(chia ~ mosaik_rpkm, main = "Chia vs Mosaik2", xlab = "Sid Mosaik", ylab = "Chia", 
    pch = 16)
legend("bottomright", "R^2=-0.06028", cex = 0.75)
abline(lm(chia ~ mosaik_rpkm))

Plot: Mosaik2 vs BWA

plot(mosaik_rpkm ~ bwa_rpkm, main = "Mosaik2 vs BWA", xlab = "Sid BWA", ylab = "Sid Mosaik", 
    pch = 16)
legend("topright", "R^2=-0.574", cex = 0.75)
abline(lm(mosaik_rpkm ~ bwa_rpkm))

Plot: Chia vs Mosaik2 Default

plot(chia ~ mosaik_def_rpkm, main = "Chia vs Mosaik2 Default", xlab = "Mosaik Default", 
    ylab = "Chia", pch = 16)
legend("bottomright", "R^2=-0.07266", cex = 0.75)
abline(lm(chia ~ mosaik_def_rpkm))

Plot: Mosaik2 Default vs BWA

plot(mosaik_def_rpkm ~ bwa_rpkm, main = "Mosaik2 Default vs BWA", xlab = "BWA", 
    ylab = "Mosaik Default", pch = 16)
legend("topright", "R^2=-0.564", cex = 0.75)
abline(lm(mosaik_rpkm ~ bwa_rpkm))

Plot: Mosaik2 Default vs Mosaik2

plot(mosaik_rpkm ~ mosaik_def_rpkm, main = "Mosaik2 Default vs Mosaik2", xlab = "Mosaik Default", 
    ylab = "Mosaik", pch = 16)
legend("bottomright", "R^2=0.9875", cex = 0.75)
abline(lm(mosaik_rpkm ~ mosaik_def_rpkm))

END