差异基因密度分布

作者: 生信云笔记 | 来源:发表于2023-03-25 12:59 被阅读0次

R语言可视化（五）：密度分布图绘制
绘制正态分布概率密度函数
多元正态分布的极大似然估计
2018-10-15
指数分布
seaborn常见绘图学习总结（回归图）
Logistic回归
dnorm pnorm qnorm rnorm详解
四、正态分布
R- 二项分布

最近做了一个似乎有些无聊的事情，就是分析RNA-seq的差异基因在染色体上的密度分布。为什么说无聊呢？因为不知道这个分析能说明什么问题，也没有看到过类似如差异基因在染色体上的密度跟什么生物意义有联系的研究。或许是咱孤陋寡闻，也或许是目前没有发现这方面的意义吧！不想那么多了，先做了再说。

计算

将染色体按照一定大小划分成bin，然后统计每个bin里面有多少个差异基因，这个计算起来不难，用R写了一个小函数来完成：

calc_density <- function(genebed, chrlen, binsize) {
  genebed <- genebed[,c(1:4)]
  colnames(genebed) <- c("chr", "start", "end", "gene")
  genebed <- genebed[order(genebed$chr), ]

  colnames(chrlen) <- c('chr','len')
  chrname <- chrlen$chr
  chrlen <- as.numeric(chrlen$len)
  names(chrlen) <- chrname

  gene_density <- list()
  for (chr in names(chrlen)) {
    num_bins <- ceiling(chrlen[chr] / binsize)
    gene_density[[chr]] <- rep(0, num_bins)
  }

  for (chr in names(chrlen)) {
    chr_genes <- genebed[genebed$chr == chr, ]

    for (i in 1:length(gene_density[[chr]])) {
      bin_start <- (i - 1) * binsize + 1
      bin_end <- min(i * binsize, chrlen[chr])
      bin_length <- bin_end - bin_start + 1

      bin_genes <- chr_genes[chr_genes$end >= bin_start & chr_genes$start <= bin_end, ]
      num_genes <- nrow(bin_genes)

      if (num_genes > 0) {
        gene_density[[chr]][i] <- num_genes
      }
    }
  }

  density_df <- data.frame()
  for (chr in names(chrlen)) {
    density_chr <- gene_density[[chr]]
    start_positions <- seq(1, chrlen[chr], by = binsize)
    end_positions <- c(start_positions[2:length(start_positions)] - 1, chrlen[chr])
    chr_df <- data.frame(Chr=chr, Start=start_positions End=end_positions, Value=density_chr)
    density_df <- rbind(density_df, chr_df)
  }

  return(density_df)
}

genebed <- read.table('dge.bed',stringsAsFactors=F,sep='\t',header=F)
head(dge)
     V1        V2        V3    V4 V5
1  chr3  90476126  90488379  Ilf2
2  chr6  28475139  28935162  Snd1
3  chr6 128362967 128376146 Foxm1
4 chr11  78301529  78322457 Spag5
5 chr17  35833921  35838306 Tubb5
6 chr17  56303321  56323486 Uhrf1

chrlen <- read.table('chrlen.txt',stringsAsFactors=F,sep='\t',header=F)
head(chrlen)
    V1        V2
1 chr1 195471971
2 chr2 182113224
3 chr3 160039680
4 chr4 156508116
5 chr5 151834684
6 chr6 149736546

genden <- calc_density(genebed, chrlen, 1000000)
head(genden)
   Chr   Start End     Value
1 chr1       1 1000000     0
2 chr1 1000001 2000000     0
3 chr1 2000001 3000000     0
4 chr1 3000001 4000000     0
5 chr1 4000001 5000000     4
6 chr1 5000001 6000000     1

该函数需要三个参数：1、genebed差异基因的位置信息的数据框，2、chrlen染色体长度数据框，3、binsize划分bin的大小。

分布

绘图使用的是RIdeogram包，ideogram函数就可以搞定了，画基因密度至少需要两个参数：1、karyotype染色体信息，至少三列不含中心粒只有染色体的起始和终止位置信息，五列的话后面两列是中心粒位置信息，2、overlaid基因密度信息，四列的数据框，包含基因位置和数据，列名分别为Chr、Start、End、Value。下面用内置的数据来演示：

library(RIdeogram)

# 内置人的染色体信息
data(human_karyotype)
head(human_karyotype)
  Chr Start       End  CE_start    CE_end
1   1     0 248956422 122026459 124932724
2   2     0 242193529  92188145  94090557
3   3     0 198295559  90772458  93655574
4   4     0 190214555  49712061  51743951
5   5     0 181538259  46485900  50059807
6   6     0 170805979  58553888  59829934

# 内置基因密度信息
data(gene_density)
head(gene_density)
  Chr   Start     End Value
1   1       1 1000000    65
2   1 1000001 2000000    76
3   1 2000001 3000000    35
4   1 3000001 4000000    30
5   1 4000001 5000000    10
6   1 5000001 6000000    10

ideogram(karyotype = human_karyotype, overlaid = gene_density)

结果如下：