genomation

Question

[daffodil]

Hi every body,

I wrote this code in the R, but I got this error:

Error in data.frame(peakID = mcols(associated_peaks)$metadata, distance = distances,  : 
  arguments imply differing number of rows: 0, 22167

It would be appreciated if you could help me.

library(Biostrings)
library(genomation)

# Assuming peak_data contains peak information
peak_data <- read.table("Desktop/ATAC data/Peak Leptonema:Zygonema Danko /SRR21230488_peak_peaks.xls", header = TRUE)

# Assuming tss_data contains TSS information
tss_data <- read.table("Desktop/Meiosis_Genes.bed", header = FALSE)

# Create GRanges objects from the data frames
pk1.gr <- GRanges(
  seqnames = Rle(peak_data$chr),
  ranges = IRanges(start = peak_data$start, end = peak_data$end),
  abs_summit = peak_data$abs_summit,  
  pileup = peak_data$pileup,          
  X.log10.pvalue. = peak_data$X.log10.pvalue.,  
  fold_enrichment = peak_data$fold_enrichment,  
  X.log10.qvalue. = peak_data$X.log10.qvalue.  
)
tss.gr <- GRanges(
  seqnames = Rle(tss_data$V1),
  ranges = IRanges(start = tss_data$V2, end = tss_data$V3)
)

unique_pk1_chr <- unique(seqnames(pk1.gr))
unique_tss_chr <- unique(seqnames(tss.gr))
# Identify common chromosomes
common_chromosomes <- intersect(unique_pk1_chr, unique_tss_chr)


# Update chromosome names in both objects using pruning.mode="coarse"
pk1.gr <- keepSeqlevels(pk1.gr, common_chromosomes, pruning.mode = "coarse")
tss.gr <- keepSeqlevels(tss.gr, common_chromosomes, pruning.mode = "coarse")

# Print the updated unique chromosome names
cat("Updated Chromosome names in pk1.gr:", unique(seqnames(pk1.gr)), "\n")
cat("Updated Chromosome names in tss.gr:", unique(seqnames(tss.gr)), "\n")

# Get the peaks that overlap with TSS regions
overlapping_peaks <- subsetByOverlaps(pk1.gr, tss.gr)

# Print or further analyze the overlapping peaks
print(overlapping_peaks)
#Calculate the count 
counts <- countOverlaps(pk1.gr, tss.gr)
head(counts)
findOverlaps(pk1.gr,tss.gr)
# find nearest CpGi to each TSS
n.ind=nearest(pk1.gr,tss.gr)
# get distance to nearest
dists=distanceToNearest(pk1.gr,tss.gr,select="arbitrary")
dists
head(dists)

# Extract distances and associated peaks
distances <- mcols(dists)$distance
associated_peaks <- peaks[queryHits(dists)]
# Check the dimensions of the 'peaks' object
dim(peaks)

# Print the indices obtained from queryHits(dists)
print(queryHits(dists))

# Create a data frame for better visualization
result_df <- data.frame(
  peakID = mcols(associated_peaks)$metadata,
  distance = distances,
  seqnames = seqnames(associated_peaks),
  start = start(associated_peaks),
  end = end(associated_peaks)
)