Merge pull request #46 from CostaLab/devel

Updated vignettes
CostaLab · Nov 12, 2024 · 032dd81 · 032dd81
2 parents e880916 + 4802126
commit 032dd81
Show file tree

Hide file tree

Showing 247 changed files with 684 additions and 1,396 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: sigurd
 Type: Package
 Title: Single cell Genotyping Using RNA Data
-Version: 0.3.6
+Version: 0.3.8
 Authors@R: c(
   person(given = "Martin",
          family = "Grasshoff",

diff --git a/NAMESPACE b/NAMESPACE
diff --git a/R/AllelFrequencyFoldChange.R b/R/AllelFrequencyFoldChange.R
diff --git a/R/AmpliconSupplementing.R b/R/AmpliconSupplementing.R
diff --git a/R/CalculateAlleleFrequency.R b/R/CalculateAlleleFrequency.R
diff --git a/R/CalculateAltReads.R b/R/CalculateAltReads.R
diff --git a/R/CalculateConsensus.R b/R/CalculateConsensus.R
diff --git a/R/CalculateCorrelationPValue.R b/R/CalculateCorrelationPValue.R
@@ -3,12 +3,13 @@
 #'@description
 #'We perform the correlation of SNVs and calculate the P values.
 #'@importFrom stats cor.test
-#'@param variant_values The fraction values you are analysing. A vector.
+#'@param variant_values The fraction values you are analyzing. A vector.
 #'@param other_mutation All other variants you have. A vector of variant names.
 #'@param all_variants_list List of fraction values for all the variants you want to compare your variant with.
 #'@param min_intersecting_cells Minimum number of intersecting cells. Correlations with less than this will not be performed.
+#'@param value_type Are we using consensus or other information?
 #'@export
-CalculateCorrelationPValue <- function(variant_values, other_mutation, all_variants_list, min_intersecting_cells = 5){
+CalculateCorrelationPValue <- function(variant_values, other_mutation, all_variants_list, value_type = "consensus", min_intersecting_cells = 5){
   other_variant_values <- all_variants_list[[other_mutation]]
   if(sum(names(variant_values) %in% names(other_variant_values)) == 0){
     #print("No intersection")
@@ -28,10 +29,17 @@ CalculateCorrelationPValue <- function(variant_values, other_mutation, all_varia
       return(result)
     } else if(length(variant_values) > 2){
       result <- stats::cor.test(variant_values, other_variant_values)
-      cells_som_alt <- sum(variant_values == 1)
-      cells_som_ref <- sum(variant_values == 0)
-      cells_MT_alt  <- sum(other_variant_values == 1)
-      cells_MT_ref  <- sum(other_variant_values == 0)
+      if(value_type == "consensus"){
+        cells_som_alt <- sum(variant_values == 1)
+        cells_som_ref <- sum(variant_values == 0)
+        cells_MT_alt  <- sum(other_variant_values == 1)
+        cells_MT_ref  <- sum(other_variant_values == 0)
+      } else{
+        cells_som_alt <- sum(variant_values > 0)
+        cells_som_ref <- sum(variant_values == 0)
+        cells_MT_alt  <- sum(other_variant_values > 0)
+        cells_MT_ref  <- sum(other_variant_values == 0)
+      }
       result <- c(result$p.value, result$estimate, cells_som_alt, cells_som_ref, cells_MT_alt, cells_MT_ref)
     } else{
       #print("We do not have more than 2 cells for the somatic variant.")

diff --git a/R/CalculateCoverage.R b/R/CalculateCoverage.R
diff --git a/R/CalculateFisherTestPValue.R b/R/CalculateFisherTestPValue.R
diff --git a/R/CalculateQuality.R b/R/CalculateQuality.R
diff --git a/R/CalculateRefReads.R b/R/CalculateRefReads.R
diff --git a/R/CalculateStrandCorrelation.R b/R/CalculateStrandCorrelation.R
diff --git a/R/CallSupport.R b/R/CallSupport.R
diff --git a/R/ClonalDefinition.R b/R/ClonalDefinition.R
diff --git a/R/ClonalDiversity.R b/R/ClonalDiversity.R
diff --git a/R/CombineSEobjects.R b/R/CombineSEobjects.R
diff --git a/R/Enrichment_Fisher.R b/R/Enrichment_Fisher.R
diff --git a/R/Filtering.R b/R/Filtering.R
diff --git a/R/GetCellInfoPerVariant.R b/R/GetCellInfoPerVariant.R
diff --git a/R/GetVariantInfo.R b/R/GetVariantInfo.R
diff --git a/R/HeatmapVOI.R b/R/HeatmapVOI.R
diff --git a/R/LoadingMAEGATK_typewise.R b/R/LoadingMAEGATK_typewise.R
diff --git a/R/LoadingRawMatrix_typewise.R b/R/LoadingRawMatrix_typewise.R
diff --git a/R/LoadingVCF_typewise.R b/R/LoadingVCF_typewise.R
diff --git a/R/LoadingVarTrix_typewise.R b/R/LoadingVarTrix_typewise.R
diff --git a/R/Merging_SE_list.R b/R/Merging_SE_list.R
diff --git a/R/RowWiseSplit.R b/R/RowWiseSplit.R
diff --git a/R/SeparatingMatrixToList.R b/R/SeparatingMatrixToList.R
diff --git a/R/SetVariantInfo.R b/R/SetVariantInfo.R
diff --git a/R/UMIQC.R b/R/UMIQC.R
diff --git a/R/UMI_seq.R b/R/UMI_seq.R
diff --git a/R/VariantBurden.R b/R/VariantBurden.R
diff --git a/R/VariantCloneSizeThresholding.R b/R/VariantCloneSizeThresholding.R
diff --git a/R/VariantCorrelationHeatmap.R b/R/VariantCorrelationHeatmap.R
diff --git a/R/VariantFisherTestHeatmap.R b/R/VariantFisherTestHeatmap.R
diff --git a/R/VariantQuantileThresholding_Combined.R b/R/VariantQuantileThresholding_Combined.R
diff --git a/R/VariantSelection_Group.R b/R/VariantSelection_Group.R
diff --git a/R/VariantSelection_Quantile.R b/R/VariantSelection_Quantile.R
@@ -22,24 +22,24 @@ VariantSelection_Quantile <- function(SE, min_coverage = 2, quantiles = c(0.1, 0
     SummarizedExperiment::assays(SE)[["fraction"]][nocall_check] <- NA
     SummarizedExperiment::assays(SE)[["coverage"]][nocall_check] <- NA
   }
-  
+
   if(verbose) print("Get the mean allele frequency and coverage.")
   mean_af  <- Matrix::rowMeans(SummarizedExperiment::assays(SE)[["fraction"]], na.rm = TRUE)
   mean_cov <- Matrix::rowMeans(SummarizedExperiment::assays(SE)[["coverage"]], na.rm = TRUE)
-  
+
   if(verbose) print("Get the quantiles of the VAFs of each variant.")
   quantiles <- lapply(quantiles, function(x) apply(SummarizedExperiment::assays(SE)[["fraction"]], 1, quantile, x, na.rm = TRUE))
-  
+
   # Apply min_quality filtering
   if(!is.null(min_quality)){
     vars <- data.frame(Mean_AF = mean_af, Mean_Cov = mean_cov, VariantQuality = SummarizedExperiment::rowData(SE)$VariantQuality, Quantile1 = quantiles[[1]], Quantile2 = quantiles[[2]])
     vars <- subset(vars, !is.na(vars$VariantQuality) & vars$VariantQuality > min_quality)
   } else{
     vars <- data.frame(Mean_AF = mean_af, Mean_Cov = mean_cov, Quantile1 = quantiles[[1]], Quantile2 = quantiles[[2]])
   }
-  
+
   if(verbose) print("Thresholding using the quantile approach.")
-  vois <- subset(vars, vars$Mean_AF > mean_allele_frequency & vars$Mean_Cov > min_coverage & vars$Quantile1 < thresholds[1] & vars$Quantile2 > thresholds[2])
-  
+  vois <- subset(vars, vars$Mean_AF > mean_allele_frequency & vars$Mean_Cov > min_coverage & vars$Quantile1 <= thresholds[1] & vars$Quantile2 >= thresholds[2])
+
   return(rownames(vois))
 }
diff --git a/R/VariantSelection_TopCells.R b/R/VariantSelection_TopCells.R
diff --git a/R/VariantSelection_VMR.R b/R/VariantSelection_VMR.R
@@ -18,7 +18,20 @@ VariantSelection_VMR <- function(SE, stabilize_variance = TRUE, low_coverage_thr
   fraction[is.na(fraction)] <- 0
   mean_af <- Matrix::rowSums(SummarizedExperiment::assays(SE)[["alts"]]) / Matrix::rowSums(SummarizedExperiment::assays(SE)[["coverage"]])
   mean_af[is.na(mean_af)] <- 0
-
+
+  if(verbose) print("We calculate the concordance.")
+  reads_forward <- SummarizedExperiment::assays(SE)[["forward"]]
+  reads_reverse <- SummarizedExperiment::assays(SE)[["reverse"]]
+  dt <- merge(data.table::data.table(summary(reads_forward)), 
+              data.table::data.table(summary(reads_reverse)), 
+              by.x = c("i", "j"), by.y = c("i", "j"), 
+              all = TRUE)
+  dt <- dt[dt[,x.x] > 0 | dt[,x.y] > 0,]
+  dt <- data.table::data.table(variant = rownames(SE)[dt[[1]]],
+                               cell_id = dt[[2]],
+                               fw = dt[[3]], rev = dt[[4]])
+  concordance <- dt[, .(cor = suppressWarnings(stats::cor(c(fw), c(rev), method = "pearson", use = "pairwise.complete"))), by = list(variant)]
+
   if(stabilize_variance){
     if(verbose) print("We stabilize the variance.")
     coverage <- SummarizedExperiment::assays(SE)[["coverage"]]
@@ -42,19 +55,17 @@ VariantSelection_VMR <- function(SE, stabilize_variance = TRUE, low_coverage_thr
 
   if(verbose) print(paste0("We check if a cells has more than ", minimum_fw_rev_reads, " reads."))
   detected <- (SummarizedExperiment::assays(SE)[["forward"]] >= minimum_fw_rev_reads) + (assays(SE)[["reverse"]] >= minimum_fw_rev_reads)
-
+  
   voi <- data.frame(
     variant = rownames(detected),
     vmr = variants_variance / (mean_af + 0.00000000001),
     vmr_log10 = log10(variants_variance / (mean_af + 0.00000000001)),
     mean = round(mean_af, 7),
     variance = round(variants_variance, 7),
     cells_detected = Matrix::rowSums(detected == 2), # We only select cells that have enough reads for both directions.
-    strand_correlation = SummarizedExperiment::rowData(SE)$Concordance,
+    strand_correlation = concordance$cor, #SummarizedExperiment::rowData(SE)$Concordance,
     mean_coverage = Matrix::rowMeans(SummarizedExperiment::assays(SE)[["coverage"]]),
     stringsAsFactors = FALSE, row.names = rownames(detected)
   )
   return(voi)
 }
-
-
diff --git a/R/VariantWiseCorrelation.R b/R/VariantWiseCorrelation.R
@@ -8,8 +8,9 @@
 #'@param n_cores Number of cores you want to use. Numeric.
 #'@param p_value_adjustment Method for P value adjustment. See p.adjust for details.
 #'@param verbose Should the function be verbose? Default = TRUE
+#'@param value_type Are we using consensus or other information?
 #'@export
-VariantWiseCorrelation <- function(variants_list, n_cores = 1, p_value_adjustment = "fdr", verbose = TRUE){
+VariantWiseCorrelation <- function(variants_list, n_cores = 1, p_value_adjustment = "fdr", value_type = "consensus", verbose = TRUE){
   # We correlate the somatic variants with each other and the MT variants.
   # Since we have tens of thousands of MT variants, we do not correlate them with each other.
   variants <- names(variants_list)
@@ -22,7 +23,7 @@ VariantWiseCorrelation <- function(variants_list, n_cores = 1, p_value_adjustmen
     variants_values_use <- variants_list[[variant_use]]
     variants_list_use <- variants_list[names(variants_list) != variant_use]
     all_variants <- names(variants_list_use)
-    results <- parallel::mclapply(X = all_variants, CalculateCorrelationPValue, variant_values = variants_values_use, all_variants_list = variants_list_use, mc.cores = n_cores)
+    results <- parallel::mclapply(X = all_variants, CalculateCorrelationPValue, variant_values = variants_values_use, all_variants_list = variants_list_use, mc.cores = n_cores, value_type = value_type)
     results <- do.call("rbind", results)
     results <- data.frame(Variant1 = variant_use, Variant2 = all_variants, P = results[,1], Corr = results[,2],
                           Cells_1_Alt = results[,3], Cells_1_Ref = results[,4], Cells_2_Alt = results[,5], Cells_2_Ref = results[,6])
@@ -35,6 +36,10 @@ VariantWiseCorrelation <- function(variants_list, n_cores = 1, p_value_adjustmen
   if(verbose) print("We remove the negative corrlated SNPs.")
   results_total <- subset(results_total, Corr > 0)
 
+  if(verbose) print("Remove duplicative results.")
+  results_check <- apply(results_total[, c("Variant1", "Variant2")], 1, function(x) paste(sort(x), collapse = "_"))
+  results_total <- results_total[!duplicated(results_check), , drop = FALSE]
+
   if(verbose) print(paste0("Adjusting P values using ", p_value_adjustment, "."))
   results_total$P_adj <- stats::p.adjust(results_total$P, method = p_value_adjustment)
   rownames(results_total) <- NULL

diff --git a/R/VariantWiseFisherTest.R b/R/VariantWiseFisherTest.R
diff --git a/R/char_to_numeric.R b/R/char_to_numeric.R
diff --git a/R/combine_NAMES.R b/R/combine_NAMES.R
diff --git a/R/combine_SparseMatrix.R b/R/combine_SparseMatrix.R
diff --git a/R/computeAFMutMatrix.R b/R/computeAFMutMatrix.R
diff --git a/R/getAltMatrix.R b/R/getAltMatrix.R
diff --git a/R/getMutMatrix.R b/R/getMutMatrix.R
diff --git a/R/getReadMatrix.R b/R/getReadMatrix.R
diff --git a/R/getRefMatrix.R b/R/getRefMatrix.R
diff --git a/R/get_consensus.R b/R/get_consensus.R
diff --git a/R/load_object.R b/R/load_object.R
diff --git a/R/save_object.R b/R/save_object.R
diff --git a/R/sdiv.R b/R/sdiv.R
diff --git a/README.md b/README.md
@@ -33,7 +33,7 @@ The mutation data was obtained from the Sanger Institute Catalogue Of Somatic Mu
 
 ```
 
-# Current Features v0.3.6
+# Current Features v0.3.8
 
 - Loading data from VarTrix and MAEGATK.
 - Transforming the data to be compatible for joint analysis.

diff --git a/docs/404.html b/docs/404.html
diff --git a/docs/ReadMe.html b/docs/ReadMe.html
diff --git a/docs/_config.yml b/docs/_config.yml
diff --git a/docs/articles/BoneMarrow_SIGURD.html b/docs/articles/BoneMarrow_SIGURD.html
diff --git a/docs/articles/BoneMarrow_SIGURD_files/figure-html/Vis_SIGURD-1.pdf b/docs/articles/BoneMarrow_SIGURD_files/figure-html/Vis_SIGURD-1.pdf
diff --git a/docs/articles/BoneMarrow_SIGURD_files/figure-html/Vis_SIGURD-1.png b/docs/articles/BoneMarrow_SIGURD_files/figure-html/Vis_SIGURD-1.png
diff --git a/docs/articles/BoneMarrow_SIGURD_files/figure-html/unnamed-chunk-3-1.png b/docs/articles/BoneMarrow_SIGURD_files/figure-html/unnamed-chunk-3-1.png
diff --git a/docs/articles/BoneMarrow_SIGURD_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/BoneMarrow_SIGURD_files/figure-html/unnamed-chunk-4-1.png