HI ,
I've got a problem I can't solve when i used LinkPeaks. I thought my error was similar with #12 and i try to solve it by adding sep = c(":", "-") but i failed.
I found some DEs and want to link the peaks and i got this error in some genes:
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
No significant links found
Testing 1 genes and 132282 peaks
Testing 1 genes and 132282 peaks
Error in .get_data_frame_col_as_numeric(df, granges_cols[["start"]]) :
some values in the "start" column cannot be turned into numeric values
Calls: CoveragePlot ... makeGRangesFromDataFrame -> .get_data_frame_col_as_numeric
In addition: Warning messages:
1: Removed 2 rows containing missing values (position_stack).
2: Removed 1 rows containing missing values (position_stack).
3: Removed 1 rows containing missing values (position_stack).
4: Removed 2 rows containing missing values (position_stack).
5: Removed 1 rows containing missing values (position_stack).
6: Expected 3 pieces. Missing pieces filled with `NA` in 1 rows [1].
Execution halted
(the error gene is HLA-DPB1)
My code:
（1）callpeak:
##call peaks using MACS2
peaks <- CallPeaks(pbmc, macs2.path = "/share/work4/wangshuang179/software/conda/Miniconda3/Miniconda3/bin/macs2")
head(peaks)
head(rownames(peaks))
peaks <- GRangesToString(peaks,sep = c(":", "-")) ##add for this error
peaks <- StringToGRanges(peaks,sep = c(":", "-")) ##add for this error
##remove peaks on nonstandard chromosomes and in genomic blacklist regions
peaks <- keepStandardChromosomes(peaks, pruning.mode = "coarse")
peaks <- subsetByOverlaps(x = peaks, ranges = blacklist_hg38_unified, invert = TRUE)
#peaks <- subsetByOverlaps(x = peaks, ranges = blacklist, invert = TRUE)
##quantify counts in each peak
peaks <- GRangesToString(peaks,sep = c(":", "-")) ##add for this error
peaks <- StringToGRanges(peaks,sep = c(":", "-")) ##add for this error
head(rownames(peaks))
#macs2_counts <- FeatureMatrix(fragments = Fragments(pbmc),features = peaks,sep = c(":", "-"),cell = colnames(pbmc),verbose = TRUE)
macs2_counts <- FeatureMatrix(fragments = Fragments(pbmc),features = peaks,sep = c(":", "-"),cell = colnames(pbmc),verbose = TRUE)
##create a new assay using the MACS2 peak set and add it to the Seurat object
pbmc[["peaks"]] <- CreateChromatinAssay(counts = macs2_counts,sep = c(":", "-"),fragments = fragpath,annotation = annotation)
**(2)LinkPeaks:**
pbmc <- readRDS(argv$rds)
##find diff gene
DefaultAssay(pbmc) <- "RNA"
markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
write.table(markers,file=paste(argv$outdir,'/',argv$sample,"_diff.xls",sep=""),sep="\t",row.names = FALSE,col.names = TRUE,quote=F)
diff <- markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
write.table(diff,file=paste(argv$outdir,'/',argv$sample,"_diff_top2.xls",sep=""),sep="\t",row.names = FALSE,col.names = TRUE,quote=F)
print("find diff gene ok")
##link peaks to genes
peak_anno_gene <- unique(as.vector(pbmc@assays$peaks@annotation$gene_name))
gene_diff <- unique(diff$gene)
GENE <- intersect(peak_anno_gene,gene_diff)
print(GENE)
cluster <- unique(pbmc@meta.data$seurat_clusters)
print(cluster)
DefaultAssay(pbmc) <- "peaks"
for (i in GENE){
print(i)
pbmc_plot <- pbmc
DefaultAssay(pbmc_plot) <- "peaks"
pbmc_plot <- RegionStats(pbmc_plot, genome = BSgenome.Hsapiens.UCSC.hg38)
head(rownames(pbmc_plot@assays$peaks))
head(rownames(pbmc_plot@assays$SCT))
pbmc_plot <- LinkPeaks(object = pbmc_plot,peak.assay = "peaks",expression.assay = "SCT",genes.use = i)
pdf(paste(argv$outdir,'/',argv$sample,"_",i,'_LinkPeaks.pdf',sep=''),onefile=F,width=6,height=6)
p1 <- CoveragePlot(object = pbmc_plot,region = i,features = i,expression.assay = "SCT",idents = cluster)
print(p1)
dev.off()
}
Some output:
(1)head(peaks):
(2)head(rownames(peaks)):
NULL
I dont known why its NULL and it looks work well with some genes.
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4
