利用R计算特定域区域的基因频率

Question

我正在研究植物的基因组拓扑学，以便了解核基因组的组织-比较化。

我想知道一种快速的方法来计算特定的基因组区域的基因密度，以便绘制该区域开始和结束的每bin +- 50Kb的基因频率。

有人有什么想法来解决这个问题？先谢谢你

Answer 1

如果你想我已经在Biostar的帖子与更多的信息，我做的。

在BioStar上发表的帖子，并附上所需输出的数字示例。

所以对于基因密度，我指的是与我的区域重叠的基因数量。但我必须将这个计数与我的区域大小进行标准化，因为它们的大小不一样。

这里你有我的区域数据框架对象。

> head(regionFileCustom)
            chr startCust  endCust                     id score strand
1 CMiso1.1chr01    -50000  3550000  Domain-1CMiso1.1chr01  1000      +
2 CMiso1.1chr01  14450000 15300000 Domain-17CMiso1.1chr01  1000      +
3 CMiso1.1chr01  15250000 15500000 Domain-18CMiso1.1chr01  1000      +
4 CMiso1.1chr01  15600000 16650000 Domain-19CMiso1.1chr01  1000      +
5 CMiso1.1chr01  16600000 17050000 Domain-20CMiso1.1chr01  1000      +
6 CMiso1.1chr01  26950000 27450000 Domain-34CMiso1.1chr01  1000      +

然后是我的基因数据框架对象。

> head(features)
            chr start   end                    id score strand
1 CMiso1.1chr01   955 24594 CMiso1.1chr01g0058061     .      +
2 CMiso1.1chr01 24595 25591 CMiso1.1chr01g0058071     .      -
3 CMiso1.1chr01 25797 30580 CMiso1.1chr01g0058081     .      +
4 CMiso1.1chr01 31006 35465 CMiso1.1chr01g0058091     .      -
5 CMiso1.1chr01 41322 41398 CMiso1.1chr01g0058101     .      -
6 CMiso1.1chr01 42694 46218 CMiso1.1chr01g0058111     .      -

所以，为了进行计数，我已经尝试了这个自己的构建函数。

# function to compute the gene density per bin and the gene frequency (gene per bin per nbr of bin per genome region)

width = 201L
step=201L
n = 100
minoverlap=21L
maxgap=-1L
ignore.strand=TRUE # I don't know yet if I have to take in count the strand for this part of the analysis.

# For test windowBin.gr = windowBinVarName

computeFeatureFreq <- function(features, regionFile, outDFName){
  # convert feature as GRange
  features.gr = makeGRangesFromDataFrame(features)
  # convert regionFile as Grange
  regionFile.gr = makeGRangesFromDataFrame(regionFile)
  # Counting gene density per genome region
  # feature.density.windowBin = countOverlaps(windowBin.gr,features.gr,ignore.strand=ignore.strand, minoverlap=minoverlap, maxgap=maxgap)
  # feature.density.windowBin = countOverlaps(features.gr, windowBin.gr,ignore.strand=ignore.strand, minoverlap=minoverlap, maxgap=maxgap)
  feature.density.region = countOverlaps(features.gr, regionFile.gr)
  # retreive coordonates for each windowBin
  regionFile.df=as.data.frame(regionFile.gr)
  feature.density.region.df = data.frame(chr=regionFile.df$seqnames, start=regionFile.df$start, end=regionFile.df$end, count=feature.density.region)
  return(feature.density.region.df)
}

如果有什么不清楚的地方，请告诉我，谢谢。