Extrating out the getMarkerFeatures

[mano2991]

trafficstars

I have two samples in my ArchR file. I need to compare them. So I used the getMarkerFeatures to do a pairwise comparison in my data. This gives me Log2Fold change, FDR, Pvalue, Mean, AUC, etc I can plot a heatmap, and Volcanotplot using these values while extracting the data. using getPeaks for the peak file and getMarkers for the genescore matrix I'm not able to extract the information it always gives me the FDR value along with the Gene Name for the genescore Matrix file and Peaks it gives me the peak locating along with FDR and Log2Fold change hence I made a custom script.

Is there any other way to extract the information instead of manually?

markersPeaks <- getMarkerFeatures(
    ArchRProj = proj5, 
    useMatrix = "PeakMatrix", 
    groupBy = "Sample",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon",
  useGroups = "Uninjured",
  bgdGroups = "Injured_7dpi"
)

A custom script for extracting the Peak file

# Load necessary library
library(SummarizedExperiment)

# Extract necessary data from elementMetadata
element_metadata <- elementMetadata(markersPeaks)
seqnames <- as.character(element_metadata$seqnames)
start <- element_metadata$start
end <- element_metadata$end

# Extract assay data
log2FC <- markersPeaks@assays@data@listData$Log2FC$x
mean <- markersPeaks@assays@data@listData$Mean$x
fdr <- markersPeaks@assays@data@listData$FDR$x
pval <- markersPeaks@assays@data@listData$Pval$x
meanDiff <- markersPeaks@assays@data@listData$MeanDiff$x
auc <- markersPeaks@assays@data@listData$AUC$x
meanBGD <- markersPeaks@assays@data@listData$MeanBGD$x

# Check the length of extracted data to ensure they are consistent
lengths <- sapply(list(seqnames, start, end, log2FC, mean, fdr, pval, meanDiff, auc, meanBGD), length)
names(lengths) <- c("seqnames", "start", "end", "log2FC", "mean", "fdr", "pval", "meanDiff", "auc", "meanBGD")
print(lengths)

# Combine extracted values into a data frame
if (all(lengths == lengths[1])) {
  peak_data <- data.frame(
    seqnames = seqnames,
    start = start,
    end = end,
    Log2FC = log2FC,
    Mean = mean,
    FDR = fdr,
    Pval = pval,
    MeanDiff = meanDiff,
    AUC = auc,
    MeanBGD = meanBGD
  )

  # Write the data frame to a CSV file
  write.csv(peak_data, "markersPeaks.csv", row.names = FALSE)
} else {
  stop("The lengths of the extracted data do not match.")
}

A custom script for the extracting the gene information.

# Extract seqnames with the correct length
seqnames <- rep(as.character(markersPeaks_GS@elementMetadata@listData$seqnames@values), markersPeaks_GS@elementMetadata@listData$seqnames@lengths)

# Combine genomic information
genomic_info <- data.frame(
  seqnames = seqnames,
  start = markersPeaks_GS@elementMetadata@listData$start,
  end = markersPeaks_GS@elementMetadata@listData$end,
  name = markersPeaks_GS@elementMetadata@listData$name
)

# Extract assay data
log2fc <- markersPeaks_GS@assays@data@listData$Log2FC$x
mean <- markersPeaks_GS@assays@data@listData$Mean$x
fdr <- markersPeaks_GS@assays@data@listData$FDR$x
pval <- markersPeaks_GS@assays@data@listData$Pval$x
mean_diff <- markersPeaks_GS@assays@data@listData$MeanDiff$x
auc <- markersPeaks_GS@assays@data@listData$AUC$x
mean_bgd <- markersPeaks_GS@assays@data@listData$MeanBGD$x

# Combine all data into a single data frame
combined_data <- data.frame(
  genomic_info,
  Log2FC = log2fc,
  Mean = mean,
  FDR = fdr,
  Pval = pval,
  MeanDiff = mean_diff,
  AUC = auc,
  MeanBGD = mean_bgd
)

# Write the combined data frame to a CSV file
write.csv(combined_data, file = "Gene_score_Pairwise_patz1_7dpi.csv", row.names = FALSE)