finished factor analysis section, started zinba wave

.gitignore

@@ -5,3 +5,4 @@ material_stat_methods/
 SRP022054/
 rse_gene.Rdata
 stockori_tmp.Rmd
+factanal.R

data/create_dds_batch_with_geno_confounding.R

@@ -39,16 +39,17 @@ dds_batch <- DESeq(dds_batch, betaPrior = TRUE)
 
 results(dds_batch)[50,]
 
+summary(results(dds_batch))
 ########################################################################
 
 ## add genotype confounding
-
-gen <- c(0,1,1,1,0,2,2,2)
+set.seed(292827)
+gen <- c(0,1,1,1,1,1,1,1,0,0,2,2,2,2,2,2)
 # gen <- 2^gen
 colData(dds_batch)$gen_back <- factor(gen)
 res <- results(DESeq(dds_batch))
 summary(res)
-idx_fc <- abs(res$log2FoldChange) > 0.8
+idx_fc <- abs(res$log2FoldChange) > 0.2
 cts <- counts(dds_batch)
 
 # adapt fold changes to reflect genetic bias 
@@ -56,20 +57,18 @@ cts <- counts(dds_batch)
 for(i in which(idx_fc)){
   
   if(res[i, "log2FoldChange"] > 0){
-    cts[i, colData(dds_batch)$condition == "B"] <- as.integer(round(cts[i,  
-           colData(dds_batch)$condition == "B"] / 2^(res[i, "log2FoldChange"]
+    cts[i, ] <- as.integer(round(cts[i, ] * 2^(res[i, "log2FoldChange"]
            / (2^gen))))
     
   } else {
-    cts[i, colData(dds_batch)$condition == "A"] <- as.integer(round(cts[i,  
-           colData(dds_batch)$condition == "A"] * 2^(res[i, "log2FoldChange"]
+    cts[i, ] <- as.integer(round(cts[i, ] * 2^(res[i, "log2FoldChange"]
            / (2^gen))))
   }
   
 }
 
 
-set.seed(292827)
+#set.seed(292827)
 # poisson <- Vectorize(function(x) {rpois(1, x)})
 # c2 <- counts(dds_batch)[idx_fc] +  t(base::apply(round(counts(dds_batch) / 1000)* gen,
 #                                        1, poisson))
@@ -79,6 +78,13 @@ dds_batch <- estimateSizeFactors(dds_batch)
 dds_batch <- DESeq(dds_batch)
 summary(results(dds_batch))
 
+# check correlations
+
+cor(as.numeric(colData(dds_batch)$condition), 
+    as.numeric(colData(dds_batch)$sex))
+
 
+cor(as.numeric(colData(dds_batch)$condition), 
+    as.numeric(colData(dds_batch)$gen_back))
 
-save(dds_batch, file = "dds_batch.RData")
\ No newline at end of file
+save(dds_batch, file = "dds_batch.RData")

factor_ana_testing_ml.R

@@ -38,6 +38,7 @@ library("pheatmap")
 library("tidyverse")
 library("Rtsne")
 library("limma")
+library("ggthemes")
 
 theme_set(theme_gray(base_size = 18))
 
@@ -222,10 +223,10 @@ sessionInfo()
 
 
 ## ----unloaAll, echo=FALSE, message=FALSE, eval = FALSE-------------------
-   
-   pkgs <- loaded_packages() %>%
-           filter(package != "devtools") %>%
-           {.$path}
-   
-   walk(pkgs, unload)
+## 
+## pkgs <- loaded_packages() %>%
+##         filter(package != "devtools") %>%
+##         {.$path}
+## 
+## walk(pkgs, unload)