最近做了一个似乎有些无聊的事情，就是分析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)
结果如下：
该包在可视化方面做的还是挺好的，还可以用来展示基因组共线性区域，网上也有专门介绍这方面的帖子，感兴趣的可以自行搜索，这里就不再赘述了。
结束语
基因密度分布用来展示差异基因在染色体上大致的位置分布，可以直观地从图上看到每个染色体区域差异基因的数量情况。但仔细一想好像缺了点什么，或许就像开头说的那样，不是很了然做基因密度分布的目的吧。
往期回顾
R绘图配色总结
saddleplot | A/B compartments
双曲线火山图一键拿捏
ChIP-seq数据质控
ChatGPT！见证AI的力量！
