Merge pull request #9 from lvclark/VCF_export

Vcf export

Author: lvclark

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: polyRAD
 Version: 1.2
-Date: 2019-11-03
+Date: 2020-05-12
 Title: Genotype Calling with Uncertainty from Sequencing Data in Polyploids and Diploids
 Authors@R: c(person("Lindsay V.", "Clark", email = "lvclark@illinois.edu",
                     role = c("aut", "cre"), 

NAMESPACE

@@ -16,9 +16,9 @@ GetLocDepth, GetLoci, GetProbableGenotypes, GetRecurrentParent,
 GetTaxa, GetWeightedMeanGenotypes, HindHe, HindHeMapping,
 InbreedingFromHindHe, IterateHWE, IterateHWE_LD,
 IteratePopStruct, IteratePopStructLD, LocusInfo, MakeTasselVcfFilter,
-MergeRareHaplotypes, MergeTaxaDepth, nTaxa, nAlleles, nLoci,
-OneAllelePerMarker, 
-PipelineMapping2Parents, RADdata, readHMC, readProcessSamMulti,
+MergeIdenticalHaplotypes, MergeRareHaplotypes, MergeTaxaDepth, nTaxa, nAlleles,
+nLoci, OneAllelePerMarker, 
+PipelineMapping2Parents, RADdata, RADdata2VCF, readHMC, readProcessSamMulti,
 readProcessIsoloci,readStacks, readTagDigger, readTASSELGBSv2,
 RemoveHighDepthLoci, RemoveMonomorphicLoci, RemoveUngenotypedLoci,
 SetBlankTaxa, SetContamRate, SetDonorParent, SetRecurrentParent,
@@ -28,7 +28,7 @@ TestOverdispersion, VCF2RADdata)
 importFrom("stats", "dist", "dmultinom", "pchisq", "lm", "na.omit", "dbinom",
            "hclust", "cutree", "as.dist", "cor", "residuals", "lm.fit", "sd",
            "median")
-importFrom("utils", "read.csv", "write.table", "write.csv")
+importFrom("utils", "read.csv", "write.table", "write.csv", "packageVersion")
 importFrom("fastmatch", "%fin%")
 importFrom("methods", "is")
 importFrom("Rcpp", "evalCpp")
@@ -68,6 +68,7 @@ S3method(GetWeightedMeanGenotypes, RADdata)
 S3method(HindHe, RADdata)
 S3method(HindHeMapping, RADdata)
 S3method(LocusInfo, RADdata)
+S3method(MergeIdenticalHaplotypes, RADdata)
 S3method(MergeRareHaplotypes, RADdata)
 S3method(MergeTaxaDepth, RADdata)
 S3method(nAlleles, RADdata)

NEWS

@@ -13,6 +13,8 @@ Duplicated Genomes" (isolocus_sorting.Rmd).  The functions readProcessSamMulti
 and readProcessIsoloci have been added for import of data from the Python
 scripts into polyRAD.
 
+The function RADdata2VCF has been added for export to VCF.
+
 A Python script has been added to help identify full tag sequences for alleles
 imported from TASSEL-GBSv2 using VCF2RADdata.
 
@@ -42,6 +44,8 @@ impact anything frin the user's perspective, but makes the code slightly
 less unweildy and allows parental genotype estimations to be used for other
 purposes.
 
+The MergeIdenticalHaplotypes function has been added.
+
 Internal functions for simulating gametes and calculating their frequencies have
 been corrected so that they still work when only one allele is being simulated.
 

R/RcppExports.R

@@ -37,6 +37,22 @@ InitHapAssign <- function(NMmat) {
     .Call('_polyRAD_InitHapAssign', PACKAGE = 'polyRAD', NMmat)
 }
 
+Hap2SNP <- function(haps, refhap, pos) {
+    .Call('_polyRAD_Hap2SNP', PACKAGE = 'polyRAD', haps, refhap, pos)
+}
+
+Hap2Hap <- function(haps, refhap, pos) {
+    .Call('_polyRAD_Hap2Hap', PACKAGE = 'polyRAD', haps, refhap, pos)
+}
+
+MakeGTstrings <- function(genotypes, ploidy) {
+    .Call('_polyRAD_MakeGTstrings', PACKAGE = 'polyRAD', genotypes, ploidy)
+}
+
+PrepVCFexport <- function(genotypes, alleles2loc, alleleDepth, alleleNucleotides, locTable, ploidy, asSNPs) {
+    .Call('_polyRAD_PrepVCFexport', PACKAGE = 'polyRAD', genotypes, alleles2loc, alleleDepth, alleleNucleotides, locTable, ploidy, asSNPs)
+}
+
 ThirdDimProd <- function(probs, ngen, ntaxa) {
     .Call('_polyRAD_ThirdDimProd', PACKAGE = 'polyRAD', probs, ngen, ntaxa)
 }

R/classes_methods.R

@@ -2274,3 +2274,41 @@ RemoveUngenotypedLoci.RADdata <- function(object, removeNonvariant = TRUE,
 
   return(object)
 }
+
+# Function to merge alleles with identical nucleotides, generally because
+# one tag was truncated with respect to the variable region.
+MergeIdenticalHaplotypes <- function(object, ...){
+  UseMethod("MergeIdenticalHaplotypes", object)
+}
+MergeIdenticalHaplotypes.RADdata <- function(object, ...){
+  if(!is.null(object$alleleFreq) || !is.null(object$depthSamplingPermutations)){
+    stop("Run MergeIdenticalHaplotypes before running any pipeline functions.")
+  }
+  
+  remal <- integer(0) # indices of alleles to remove
+  for(L in 1:nLoci(object)){
+    thesecol <- which(object$alleles2loc == L)
+    dup <- duplicated(object$alleleNucleotides[thesecol])
+    for(al in thesecol[dup]){
+      # find allele to merge this one into
+      alM <- min(thesecol[object$alleleNucleotides[thesecol] == 
+                            object$alleleNucleotides[al]])
+      stopifnot(al != alM)
+      # consolidate read depth
+      object$alleleDepth[,alM] <- object$alleleDepth[,alM] +
+        object$alleleDepth[,al]
+      object$antiAlleleDepth[,alM] <- object$antiAlleleDepth[,alM] -
+        object$alleleDepth[,al]
+    }
+    remal <- c(remal, thesecol[dup])
+  }
+  
+  # remove duplicated alleles from all slots
+  object$alleleDepth <- object$alleleDepth[,-remal]
+  object$antiAlleleDepth <- object$antiAlleleDepth[,-remal]
+  object$alleles2loc <- object$alleles2loc[-remal]
+  object$alleleNucleotides <- object$alleleNucleotides[-remal]
+  object$depthRatio <- object$depthRatio[,-remal]
+  
+  return(object)
+}

R/data_export.R

@@ -315,3 +315,107 @@ Export_GWASpoly <- function(object, file, naIfZeroReads = TRUE){
   # export
   write.csv(outdata, file = file, row.names = FALSE, quote = FALSE)
 }
+
+RADdata2VCF <- function(object, file = NULL, asSNPs = TRUE, hindhe = TRUE,
+                        sampleinfo = data.frame(row.names = GetTaxa(object)),
+                        contigs = data.frame(row.names = unique(object$locTable$Chr))){
+  # shortcuts to functions to use
+  DataFrame <- S4Vectors::DataFrame
+  
+  varok <- attr(object$alleleNucleotides, "Variable_sites_only")
+  if(is.null(varok) || varok){
+    stop("Complete haplotype information not provided; unable to determine SNP positions.  Use refgenome argument in VCF2RADdata.")
+  }
+  if(is.null(object$locTable$Chr) || is.null(object$locTable$Pos)){
+    stop("Need chromosome and position information (Chr and Pos in locTable).")
+  }
+  if(is.null(object$locTable$Ref)){
+    warning("Reference allele not indicated.  Using the major allele for each locus.")
+    object$locTable$Ref <- object$alleleNucleotides[OneAllelePerMarker(object, commonAllele = TRUE)]
+  }
+  if(nrow(sampleinfo) != nTaxa(object) || !all(rownames(sampleinfo) %in% GetTaxa(object))){
+    "sampleinfo doesn't match taxa in RADdata object."
+  }
+  
+  # Determine most probable genotypes, and their ploidies
+  temp <- GetProbableGenotypes(object, omit1allelePerLocus = FALSE)
+  geno <- temp$genotypes
+  pld_index <- temp$ploidy_index
+  pld_ind_per_loc <- 
+    tapply(pld_index, object$alleles2loc,
+           function(x){
+             u <- unique(x)
+             if(length(u) == 1){
+               return(u)
+             } else {
+               return(as.integer(names(which.max(table(x)))))
+             }
+           })
+  pld_per_loc <- sapply(object$priorProbPloidies, sum)[pld_ind_per_loc]
+  
+  # Process data with internal RCpp function
+  temp <- PrepVCFexport(geno, object$alleles2loc, object$alleleDepth,
+                        object$alleleNucleotides, object$locTable, pld_per_loc,
+                        asSNPs)
+  REF <- Biostrings::DNAStringSet(temp$REF)
+  ALT <- Biostrings::DNAStringSetList(temp$ALT)
+  CHROM <- as.character(object$locTable$Chr)[temp$Lookup]
+  rr <- GenomicRanges::GRanges(CHROM,
+                IRanges::IRanges(start = temp$POS, 
+                                 width = BiocGenerics::width(REF)))
+  fixed <- DataFrame(REF = REF, ALT = ALT)
+  cd <- DataFrame(row.names = GetTaxa(object))
+  if(ncol(sampleinfo) > 0){
+    cd <- cbind(cd, sampleinfo[GetTaxa(object),])
+    colnames(cd) <- colnames(sampleinfo)
+  }
+  DP <- t(object$locDepth[,as.character(temp$Lookup)])
+  rownames(DP) <- NULL
+  info <- DataFrame(NS = rowSums(DP > 0), DP = rowSums(DP), LU = temp$Lookup)
+  infohdr <- DataFrame(row.names = c("NS", "DP", "LU"), Number = c("1", "1", "1"),
+                       Type = c("Integer", "Integer", "Integer"),
+                       Description = c("Number of samples with data", "Combined depth across samples",
+                                       "Lookup index of marker in RADdata object"))
+  metahdr <- DataFrame()
+  if(ncol(contigs) > 0){
+    ctg <- DataFrame(row.names = rownames(contigs), contigs)
+    colnames(ctg) <- colnames(contigs)
+  } else {
+    ctg <- DataFrame(row.names = rownames(contigs))
+  }
+  
+  # Add Hind/He if desired
+  if(hindhe){
+    infohdr <- rbind(infohdr,
+                     DataFrame(row.names = "HH", Number = "1", Type = "Float",
+                               Description = "Hind/He for the locus in the RADdata object"))
+    metahdr <- DataFrame(row.names = "HH", Number = "1", Type = "Float",
+                         Description = "Hind/He for the sample, averaged across loci in the RADdata object")
+    hh <- HindHe(object)
+    info$HH <- colMeans(hh, na.rm = TRUE)[temp$Lookup]
+    cd$HH <- rowMeans(hh, na.rm = TRUE)
+  }
+  
+  # Build VCF object
+  hdr <- VariantAnnotation::VCFHeader(reference = rownames(ctg),
+    samples = GetTaxa(object),
+    IRanges::DataFrameList(fileformat = DataFrame(row.names = "fileformat", Value = "VCFv4.3"),
+                           fileDate = DataFrame(row.names = "fileDate", Value = gsub("-", "", Sys.Date())),
+                           source = DataFrame(row.names = "source", Value = paste0("polyRADv", packageVersion("polyRAD"))),
+                           FORMAT = DataFrame(row.names = c("GT", "AD", "DP"),
+                                              Number = c("1", "R", "1"), Type = c("String", "Integer", "Integer"),
+                                              Description = c("Genotype", "Read depth for each allele", "Read depth")),
+                           INFO = infohdr, META = metahdr, contig = ctg))
+  if(ncol(cd) > 0){
+    VariantAnnotation::meta(hdr)$SAMPLE <- cd
+  }
+  vcf <- VariantAnnotation::VCF(rowRanges = rr, fixed = fixed, info = info, colData = cd,
+                                geno = S4Vectors::SimpleList(GT = temp$GT, AD = temp$AD, DP = DP),
+                                exptData = list(header = hdr), collapsed = TRUE)
+  
+  # output
+  if(!is.null(file)){
+    VariantAnnotation::writeVcf(vcf, file)
+  }
+  return(vcf)
+}