Merge pull request #24 from CostaLab/devel

Merge after adding Function Calling and Tests.
CostaLab · Nov 10, 2023 · 6ff092a · 6ff092a
2 parents 9ebfab3 + 3b84d77
commit 6ff092a
Show file tree

Hide file tree

Showing 189 changed files with 3,481 additions and 1,087 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: sigurd
 Type: Package
 Title: Single cell Genotyping Using RNA Data
-Version: 0.2.23
+Version: 0.2.32
 Authors@R: c(
   person(given = "Martin",
          family = "Grasshoff",
@@ -20,21 +20,38 @@ Encoding: UTF-8
 LazyData: true
 Imports: 
     archive,
-    bigmemory,
-    circlize,
+    BiocGenerics,
     ComplexHeatmap,
+    circlize,
+    data.table,
     dplyr,
     fastmatch,
+    GenomeInfoDb,
     ggplot2,
+    ggsci,
+    glue,
     grid,
-    magrittr,
     Matrix,
+    MatrixGenerics,
+    magrittr,
+    methods,
     parallel,
     rcompanion,
     S4Vectors,
     Seurat,
     SummarizedExperiment,
+    scales,
+    tibble,
     tidyr,
     tidyverse,
     VariantAnnotation
 RoxygenNote: 7.2.3
+Suggests:
+    GenomicRanges,
+    IRanges,
+    knitr,
+    rmarkdown,
+    SeuratObject,
+    testthat (>= 3.0.0)
+VignetteBuilder: knitr
+Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
@@ -15,6 +15,7 @@ export(GetCellInfoPerVariant)
 export(GetVariantInfo)
 export(HeatmapVoi)
 export(LoadingMAEGATK_typewise)
+export(LoadingVCF_typewise)
 export(LoadingVarTrix_typewise)
 export(Merging_SE_list)
 export(RowWiseSplit)
@@ -27,10 +28,12 @@ export(VariantFisherTestHeatmap)
 export(VariantQuantileThresholding)
 export(VariantWiseCorrelation)
 export(VariantWiseFisherTest)
+export(char_to_numeric)
 export(combine_NAMES)
 export(combine_SparseMatrix)
 export(computeAFMutMatrix)
 export(getAltMatrix)
+export(getMutMatrix)
 export(getReadMatrix)
 export(getRefMatrix)
 export(get_consensus)
@@ -39,26 +42,54 @@ export(load_object)
 export(save_object)
 export(sdiv)
 import(BiocGenerics)
-import(ComplexHeatmap)
-import(Matrix)
-import(MatrixGenerics)
-import(Seurat)
-import(SummarizedExperiment)
-import(VariantAnnotation)
 import(archive)
-import(assertthat)
-import(circlize)
-import(data.table)
-import(dplyr)
-import(fastmatch)
-import(ggplot2)
 import(ggsci)
-import(glue)
-import(grid)
-import(parallel)
-import(rcompanion)
-import(scales)
-import(stats)
-import(tibble)
-import(tidyr)
-import(tidyverse)
+importFrom(BiocGenerics,start)
+importFrom(ComplexHeatmap,Heatmap)
+importFrom(ComplexHeatmap,columnAnnotation)
+importFrom(ComplexHeatmap,draw)
+importFrom(ComplexHeatmap,rowAnnotation)
+importFrom(GenomeInfoDb,seqnames)
+importFrom(Matrix,colMeans)
+importFrom(Matrix,colSums)
+importFrom(Matrix,readMM)
+importFrom(Matrix,rowMeans)
+importFrom(Matrix,rowSums)
+importFrom(Matrix,sparseMatrix)
+importFrom(Matrix,summary)
+importFrom(Matrix,t)
+importFrom(S4Vectors,merge)
+importFrom(Seurat,AddMetaData)
+importFrom(SummarizedExperiment,SummarizedExperiment)
+importFrom(SummarizedExperiment,assays)
+importFrom(SummarizedExperiment,colData)
+importFrom(SummarizedExperiment,rowData)
+importFrom(SummarizedExperiment,rowRanges)
+importFrom(VariantAnnotation,alt)
+importFrom(VariantAnnotation,info)
+importFrom(VariantAnnotation,readGeno)
+importFrom(VariantAnnotation,readVcf)
+importFrom(VariantAnnotation,ref)
+importFrom(circlize,colorRamp2)
+importFrom(data.table,data.table)
+importFrom(dplyr,"%>%")
+importFrom(dplyr,left_join)
+importFrom(dplyr,sym)
+importFrom(glue,glue)
+importFrom(grDevices,dev.off)
+importFrom(grDevices,png)
+importFrom(grid,gpar)
+importFrom(grid,unit)
+importFrom(methods,as)
+importFrom(parallel,mclapply)
+importFrom(scales,hue_pal)
+importFrom(stats,cor.test)
+importFrom(stats,fisher.test)
+importFrom(stats,na.omit)
+importFrom(stats,p.adjust)
+importFrom(stats,quantile)
+importFrom(tibble,as_tibble)
+importFrom(tibble,tibble)
+importFrom(utils,read.csv)
+importFrom(utils,read.table)
+importFrom(utils,tail)
diff --git a/R/AmpliconSupplementing.R b/R/AmpliconSupplementing.R
@@ -2,83 +2,89 @@
 #'
 #'@description
 #'We replace the values from an scRNAseq experiment with values we have from an amplicon experiment.
-#'
-#'@import archive Matrix SummarizedExperiment VariantAnnotation
+#'@importFrom S4Vectors merge
+#'@importFrom SummarizedExperiment colData rowData assays SummarizedExperiment
+#'@importFrom stats na.omit
 #'@param scRNAseq The SummarizedExperiment object containing the scRNAseq data. 
 #'@param amplicon The SummarizedExperiment object containing the amplicon data.
+#'@param verbose Should the function be verbose? Default = TRUE
 #'@export
-AmpliconSupplementing <- function(scRNAseq, amplicon){
+AmpliconSupplementing <- function(scRNAseq, amplicon, verbose = TRUE){
   # We supplement the scRNAseq data with the amplicon data.
-  print("We get the new meta data.")
-  new_meta_data <- merge(colData(scRNAseq), colData(amplicon), by = "Cell", all.x = TRUE, all.y = TRUE,
-                         suffixes = c("scRNAseq", "Amplicon"))
+  if(verbose) print("We get the new meta data.")
+  new_meta_data <- S4Vectors::merge(SummarizedExperiment::colData(scRNAseq), SummarizedExperiment::colData(amplicon), by = "Cell", all.x = TRUE, all.y = TRUE,
+                                    suffixes = c("scRNAseq", "Amplicon"))
   rownames(new_meta_data) <- new_meta_data$Cell
   # We add an AverageCoverage column to the new meta data.
   new_meta_data$AverageCoverage                           <- new_meta_data$AverageCoveragescRNAseq
   amplicon_value                                          <- new_meta_data$AverageCoverageAmplicon
   names(amplicon_value)                                   <- colnames(amplicon)
-  amplicon_value                                          <- na.omit(amplicon_value)
+  amplicon_value                                          <- stats::na.omit(amplicon_value)
   new_meta_data[names(amplicon_value), "AverageCoverage"] <- amplicon_value
 
-  new_row_data <- merge(rowData(scRNAseq), rowData(amplicon), by = "VariantName", all.x = TRUE, all.y = TRUE,
+  new_row_data <- merge(SummarizedExperiment::rowData(scRNAseq), SummarizedExperiment::rowData(amplicon), by = "VariantName", all.x = TRUE, all.y = TRUE,
                         suffixes = c("scRNAseq", "Amplicon"))
   rownames(new_row_data) <- new_row_data$VariantName
-  # We add a VariantQuality column to the row data, showing the scRNAseq quality with the supplemented amplicon quality.
-  # We do the same for the concordance and the depth.
-  new_row_data$VariantQuality                           <- new_row_data$VariantQualityscRNAseq
-  amplicon_value                                        <- rowData(amplicon)$VariantQuality
-  names(amplicon_value)                                 <- rownames(amplicon)
-  amplicon_value                                        <- na.omit(amplicon_value)
-  new_row_data[names(amplicon_value), "VariantQuality"] <- amplicon_value
-
-  new_row_data$Concordance                           <- new_row_data$ConcordancescRNAseq
-  amplicon_value                                     <- rowData(amplicon)$Concordance
-  names(amplicon_value)                              <- rownames(amplicon)
-  amplicon_value                                     <- na.omit(amplicon_value)
-  new_row_data[names(amplicon_value), "Concordance"] <- amplicon_value
-
+  if("VariantQualityscRNAseq" %in% colnames(new_row_data)){
+    # We add a VariantQuality column to the row data, showing the scRNAseq quality with the supplemented amplicon quality.
+    # We do the same for the concordance and the depth.
+    new_row_data$VariantQuality                           <- new_row_data$VariantQualityscRNAseq
+    amplicon_value                                        <- SummarizedExperiment::rowData(amplicon)$VariantQuality
+    names(amplicon_value)                                 <- rownames(amplicon)
+    amplicon_value                                        <- stats::na.omit(amplicon_value)
+    new_row_data[names(amplicon_value), "VariantQuality"] <- amplicon_value
+  }
+
+  if("ConcordancescRNAseq" %in% colnames(new_row_data)){
+    new_row_data$Concordance                           <- new_row_data$ConcordancescRNAseq
+    amplicon_value                                     <- SummarizedExperiment::rowData(amplicon)$Concordance
+    names(amplicon_value)                              <- rownames(amplicon)
+    amplicon_value                                     <- stats::na.omit(amplicon_value)
+    new_row_data[names(amplicon_value), "Concordance"] <- amplicon_value
+  }
+
   new_row_data$Depth                           <- new_row_data$DepthscRNAseq
-  amplicon_value                               <- rowData(amplicon)$Depth
+  amplicon_value                               <- SummarizedExperiment::rowData(amplicon)$Depth
   names(amplicon_value)                        <- rownames(amplicon)
-  amplicon_value                               <- na.omit(amplicon_value)
+  amplicon_value                               <- stats::na.omit(amplicon_value)
   new_row_data[names(amplicon_value), "Depth"] <- amplicon_value
 
-  print("We get all cells and variants.")
+  if(verbose) print("We get all cells and variants.")
   all_cells     <- unique(c(colnames(scRNAseq), colnames(amplicon)))
   all_variants  <- unique(c(rownames(scRNAseq), rownames(amplicon)))
   new_meta_data <- new_meta_data[all_cells,]
   new_row_data  <- new_row_data[all_variants,]
 
-  print("We generate our output matrices.")
+  if(verbose) print("We generate our output matrices.")
   consensus <- matrix(0, ncol = length(all_cells), nrow = length(all_variants))
   rownames(consensus) <- all_variants
   colnames(consensus) <- all_cells
-  #consensus <- sparseMatrix(i = 1, j = 1, dims = c(length(all_variants), length(all_cells)), repr = "C")
+  #consensus <- Matrix::sparseMatrix(i = 1, j = 1, dims = c(length(all_variants), length(all_cells)), repr = "C")
   fraction <- consensus
   reads <- consensus
   alts <- consensus
   refs <- consensus
 
-  print("We fill the output matrices.")
-  consensus[rownames(scRNAseq), colnames(scRNAseq)] <- as.matrix(assays(scRNAseq)$consensus)
-  fraction[rownames(scRNAseq), colnames(scRNAseq)] <- as.matrix(assays(scRNAseq)$fraction)
-  reads[rownames(scRNAseq), colnames(scRNAseq)] <- as.matrix(assays(scRNAseq)$coverage)
-  alts[rownames(scRNAseq), colnames(scRNAseq)] <- as.matrix(assays(scRNAseq)$alts)
-  refs[rownames(scRNAseq), colnames(scRNAseq)] <- as.matrix(assays(scRNAseq)$refs)
+  if(verbose) print("We fill the output matrices.")
+  consensus[rownames(scRNAseq), colnames(scRNAseq)] <- as.matrix(SummarizedExperiment::assays(scRNAseq)$consensus)
+  fraction[rownames(scRNAseq), colnames(scRNAseq)]  <- as.matrix(SummarizedExperiment::assays(scRNAseq)$fraction)
+  reads[rownames(scRNAseq), colnames(scRNAseq)]     <- as.matrix(SummarizedExperiment::assays(scRNAseq)$coverage)
+  alts[rownames(scRNAseq), colnames(scRNAseq)]      <- as.matrix(SummarizedExperiment::assays(scRNAseq)$alts)
+  refs[rownames(scRNAseq), colnames(scRNAseq)]      <- as.matrix(SummarizedExperiment::assays(scRNAseq)$refs)
 
-  print("We add the the amplicon information.")
-  consensus[rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$consensus)
-  fraction[rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$fraction)
-  reads[rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$coverage)
-  alts[rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$alts)
-  refs[rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$refs)
+  if(verbose) print("We add the the amplicon information.")
+  consensus[rownames(amplicon), colnames(amplicon)] <- as.matrix(SummarizedExperiment::assays(amplicon)$consensus)
+  fraction[rownames(amplicon), colnames(amplicon)]  <- as.matrix(SummarizedExperiment::assays(amplicon)$fraction)
+  reads[rownames(amplicon), colnames(amplicon)]     <- as.matrix(SummarizedExperiment::assays(amplicon)$coverage)
+  alts[rownames(amplicon), colnames(amplicon)]      <- as.matrix(SummarizedExperiment::assays(amplicon)$alts)
+  refs[rownames(amplicon), colnames(amplicon)]      <- as.matrix(SummarizedExperiment::assays(amplicon)$refs)
 
-  #print("We add the the amplicon information.")
-  #assays(scRNAseq)[["consensus"]][rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$consensus)
-  #assays(scRNAseq)[["fraction"]][rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$fraction)
-  #assays(scRNAseq)[["coverage"]][rownames(amplicon), colnames(amplicon)] <- as.matrix(assays(amplicon)$coverage)
+  #if(verbose) print("We add the the amplicon information.")
+  #SummarizedExperiment::assays(scRNAseq)[["consensus"]][rownames(amplicon), colnames(amplicon)] <- as.matrix(SummarizedExperiment::assays(amplicon)$consensus)
+  #SummarizedExperiment::assays(scRNAseq)[["fraction"]][rownames(amplicon), colnames(amplicon)]  <- as.matrix(SummarizedExperiment::assays(amplicon)$fraction)
+  #SummarizedExperiment::assays(scRNAseq)[["coverage"]][rownames(amplicon), colnames(amplicon)]  <- as.matrix(SummarizedExperiment::assays(amplicon)$coverage)
 
-  se <- SummarizedExperiment(assays = list(consensus = as(consensus, "dgCMatrix"), fraction = as(fraction, "dgCMatrix"), coverage = as(reads, "dgCMatrix"), alts = as(alts, "dgCMatrix"), refs = as(refs, "dgCMatrix")),
-                             colData = new_meta_data, rowData = new_row_data)
+  se <- SummarizedExperiment::SummarizedExperiment(assays = list(consensus = methods::as(consensus, "dgCMatrix"), fraction = methods::as(fraction, "dgCMatrix"), coverage = methods::as(reads, "dgCMatrix"), alts = methods::as(alts, "dgCMatrix"), refs = methods::as(refs, "dgCMatrix")),
+                                                   colData = new_meta_data, rowData = new_row_data)
   return(se)
 }
diff --git a/R/CalculateAlleleFrequency.R b/R/CalculateAlleleFrequency.R
@@ -1,11 +1,13 @@
 #'Calculating the Minor Allele Frequency.
 #'@description
-#'We calculate the MAF for the MAEGATK results.
-#'@import MatrixGenerics SummarizedExperiment
+#'We calculate the MAF from a reference reads matrix and an alternative reads matrix.
+#'This function is intended to be used with the mitochondrial genome and not with other somatic mutations.
+# #'@import MatrixGenerics SummarizedExperiment
 #'@param reference_reads Reference reads matrix.
 #'@param alternative_reads List of matrices for the alternative reads.
+#'@param pseudo_count = What is the pseudo count you want to add to the reference_reads matrix. Default = 0
 #'@export
-CalculateAlleleFrequency <- function(reference_reads, alternative_reads){
+CalculateAlleleFrequency <- function(reference_reads, alternative_reads, pseudo_count = 0){
   # We remove the potential N at position 3107 of the human genome.
   alternative_reads <- alternative_reads[!grepl("_N", rownames(alternative_reads)),]
   # We get the first part of the ref row name. This includes the position and the ref allele.
@@ -19,7 +21,7 @@ CalculateAlleleFrequency <- function(reference_reads, alternative_reads){
   # We match the new ref reads matrices to be the same order as the alt read matrices.
   reference_reads <- reference_reads[match(rows_alt_reads, rows_ref_reads),]
   # We divide the alt matrix by alt + ref matrix.
-  allelefrequency <- as.matrix(alternative_reads / (alternative_reads + reference_reads))
+  allelefrequency <- as.matrix(alternative_reads / (alternative_reads + reference_reads + pseudo_count))
   allelefrequency[is.na(allelefrequency)] <- 0
   rownames(allelefrequency) <- gsub(">", "_", rownames(allelefrequency))
   return(allelefrequency)

diff --git a/R/CalculateAltReads.R b/R/CalculateAltReads.R
@@ -1,25 +1,26 @@
 #'CalculateAltReads
 #'@description
 #'We calculate the number of reads covering a variant using forward and reverse reads.
-#'@import MatrixGenerics SummarizedExperiment
+# #'@import MatrixGenerics
+#'@importFrom SummarizedExperiment SummarizedExperiment assays rowRanges
 #'@param SE SummarizedExperiment object.
 #'@param chromosome_prefix List of matrices for the alternative reads.
 #'@export
 CalculateAltReads <- function(SE, chromosome_prefix = "chrM"){
-  ref_allele <- as.character(rowRanges(SE)$refAllele)
-  reads_A <- assays(SE)[["A_counts_fw"]] + assays(SE)[["A_counts_rev"]]
+  ref_allele <- as.character(SummarizedExperiment::rowRanges(SE)$refAllele)
+  reads_A <- SummarizedExperiment::assays(SE)[["A_counts_fw"]] + SummarizedExperiment::assays(SE)[["A_counts_rev"]]
   rownames(reads_A) <- paste0(chromosome_prefix, "_", 1:nrow(reads_A), "_", ref_allele, "_A")
   reads_A <- reads_A[ref_allele != "A",]
 
-  reads_C <- assays(SE)[["C_counts_fw"]] + assays(SE)[["C_counts_rev"]]
+  reads_C <- SummarizedExperiment::assays(SE)[["C_counts_fw"]] + SummarizedExperiment::assays(SE)[["C_counts_rev"]]
   rownames(reads_C) <- paste0(chromosome_prefix, "_", 1:nrow(reads_C), "_", ref_allele, "_C")
   reads_C <- reads_C[ref_allele != "C",]
 
-  reads_G <- assays(SE)[["G_counts_fw"]] + assays(SE)[["G_counts_rev"]]
+  reads_G <- SummarizedExperiment::assays(SE)[["G_counts_fw"]] + SummarizedExperiment::assays(SE)[["G_counts_rev"]]
   rownames(reads_G) <- paste0(chromosome_prefix, "_", 1:nrow(reads_G), "_", ref_allele, "_G")
   reads_G <- reads_G[ref_allele != "G",]
 
-  reads_T <- assays(SE)[["T_counts_fw"]] + assays(SE)[["T_counts_rev"]]
+  reads_T <- SummarizedExperiment::assays(SE)[["T_counts_fw"]] + SummarizedExperiment::assays(SE)[["T_counts_rev"]]
   rownames(reads_T) <- paste0(chromosome_prefix, "_", 1:nrow(reads_T), "_", ref_allele, "_T")
   reads_T <- reads_T[ref_allele != "T",]