按定义的间隔中的 bp 位置值对行进行分组或分箱，并对它们的计数求和在 R 或 python 中

编程问答 2022-05-03

问题描述

我有两个小鼠品系（B6 和 CAST）的 SNP 文件，格式如下：

Chr     Position  B6.SNP CAST.SNP  B6.Count  CAST.Count
chr17   1000      A      G         102       82
chr17   2000      A      G         91        76
chr17   5000      C      T         55        87
chr17   7500      G      A         70        36
chr17   10200     A      G         83        45
chr17   17000     C      T         34        91
chr17   20000     G      A         95        46

我想将 SNP（即行）按 'Chr' 和 'Position' 分组到 10,000 bp 的染色体上；换句话说，将所有落在 0-10,000bp、10,0001-20,000bp 范围内的 SNP 分组，依此类推）。此外，对于落入每个区间内的所有 SNP，我想将它们的 B6 和 CAST 计数相加，并创建一个包含落入每个区间内的 SNP 数量的新列。

所以也许是一个输出文件（使用上面的例子）：

 Chr    Start   End     SNP.Count   B6.Count.Sum    CAST.Count.Sum
 chr17  0       10000   4           318             281
 chr17  10001   20000   3           212             182

提前致谢。

解决方法

这是在 R 中使用 GenomicRanges、TxDb 和 plyranges 的开始：

## Read snp data
snps <- read.table(text = "
             Chr     Position  B6.SNP CAST.SNP  B6.Count  CAST.Count
  chr17   1000      A      G         102       82
  chr17   2000      A      G         91        76
  chr17   5000      C      T         55        87
  chr17   7500      G      A         70        36
  chr17   10200     A      G         83        45
  chr17   17000     C      T         34        91
  chr17   20000     G      A         95        46
  
             ",header = T)

## Convert to GRanges object
library(GenomicRanges)
gr <- GRanges(seqnames = snps$Chr,ranges = IRanges(start = snps$Position,end = snps$Position))
mcols(gr) <- snps[,3:6]


## Use txdb to assign seqinfo
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene

names_txdb <- seqnames(seqinfo(txdb))
names_gr   <- seqnames(seqinfo(gr))

seqinfo(gr) <- seqinfo(txdb)[names_txdb[names_txdb %in% names_gr]]

## Make windows
windows <- tileGenome(seqinfo(gr),tilewidth = 10000,cut.last.tile.in.chrom = TRUE)

## Use plyranges to summarize by group overlap
library(plyranges)

df <- 
  windows %>%
  group_by_overlaps(gr) %>%
  summarise(B6.Count.Sum = sum(B6.Count),CAST.Count.Sum = sum(CAST.Count),SNP.Count = n())

binnedGR <- windows[df$query] %>% `mcols<-`(value = df[-1])

## Result as GRanges or as data.frame
binnedGR
as.data.frame(binnedGR)
binningdataframepythonrrsumsum