Can markers be calculated independently? about ctlmapping HOT 4 CLOSED

library(devtools)
install_github("DannyArends/CTLmapping", subdir = "Rctl")

library(qtl)
library(MASS)
library(ctl)

# Load Arabidopsis Thaliana data set
data(ath.metabolites)

# Define genotype and phenotype matrices
matrix_geno = ath.metab$genotypes
matrix_pheno = ath.metab$phenotypes

# Function mat_split from stackoverflow
# http://stackoverflow.com/questions/24299171/function-to-split-a-matrix-into-sub-matrices-in-r
mat_split <- function(M, r, c){
  nr <- ceiling(nrow(M)/r)
  nc <- ceiling(ncol(M)/c)
  newM <- matrix(NA, nr*r, nc*c)
  newM[1:nrow(M), 1:ncol(M)] <- M
  div_k <- kronecker(matrix(seq_len(nr*nc), nr, byrow = TRUE), matrix(1, r, c))
  matlist <- split(newM, div_k)
  N <- length(matlist)
  mats <- unlist(matlist)
  dim(mats)<-c(r, c, N)
  return(mats)
}

# Select the first marker and perform ctlscan
matrix_geno_one_marker_1 = as.matrix(mat_split(matrix_geno, nrow(matrix_geno), 1)[,,1])
ctlscan_one_marker_1 = CTLscan(genotypes = matrix_geno_one_marker_1 , phenotypes = matrix_pheno, adjust=FALSE)

# Select the second marker and perform ctlscan
matrix_geno_one_marker_2 = as.matrix(mat_split(matrix_geno, nrow(matrix_geno), 1)[,,2])
ctlscan_one_marker_2 = CTLscan(genotypes = matrix_geno_one_marker_2 , phenotypes = matrix_pheno, adjust=FALSE)

# Select the first and second markers and perform ctlscan
matrix_geno_two_markers_1_and_2 = as.matrix(mat_split(matrix_geno, nrow(matrix_geno), 2)[,,1])
ctlscan_two_markers_1_and_2 = CTLscan(genotypes = matrix_geno_two_markers_1_and_2, phenotypes = matrix_pheno, adjust=FALSE)

# Print phenotypes where the ctlscan results are different for the same marker but in different sets.
for(phenotype_i in 1:(ncol(matrix_pheno))){
  if (ctlscan_one_marker_1[[phenotype_i]]$ctl[1,] != ctlscan_two_markers_1_and_2[[phenotype_i]]$ctl[1,] && any(ctlscan_one_marker_1[[phenotype_i]]$ctl[1,] > 0.0)){
    print("ctlscan of marker 1 is not equivalent for scanning with one vs two markers for:")
    print(colnames(matrix_pheno)[phenotype_i])
    #print(ctlscan_one_marker_1[[phenotype_i]]$ctl[1,])
    #print(ctlscan_two_markers_1_and_2[[phenotype_i]]$ctl[1,])
  }
  if (ctlscan_one_marker_2[[phenotype_i]]$ctl[1,] != ctlscan_two_markers_1_and_2[[phenotype_i]]$ctl[2,] && any(ctlscan_one_marker_2[[phenotype_i]]$ctl[1,] > 0.0)){
    print("ctlscan of marker 2 is not equivalent for scanning with one vs two markers for:")
    print(colnames(matrix_pheno)[phenotype_i])
    #print(ctlscan_one_marker_2[[phenotype_i]]$ctl[1,])
    #print(ctlscan_two_markers_1_and_2[[phenotype_i]]$ctl[2,])
  }
}