用R画出染色体修饰图谱--超详细版本

用R画出染色体修饰图谱—染色体之间的物理长度之比不变

准备数据

1.染色体长度信息文件

全基因组文件建立index时可获取

2.将染色体以1Mb 滑动划分区间

bedtools makewindows -g genome.txt -w 1000000 > genome_windows_1Mb.bed

计算深度

1.准备mapping之后的bam文件

2.按照染色体位置进行排序

samtools sort wheat.bam > wheat.sorted.bam

3.排序后的bam文件建立Index

samtools index -c -@ 10 wheat.sorted.bam > wheat.sorted.bam.csi
ps:大数据的基因组比如小麦需选择生成“csi”格式的index文件.

4.将bam文件放到划分好的基因组区域计算深度

bedtools coverage -a genome_windows_1Mb.bed -b wheat.sorted.bam -sorted > wheat_depth.txt
ps:参数-sorted 适用于大基因组物种

表头分别为：染色体、区间起始位点、区间结束位点、区间内的reads数、区间内的碱基数、区间大小、区间的平均覆盖度

画图

// ###成功版本 get
setwd("D:/A学习库/小麦/reads depth in chr/")
library(ggplot2)
data <- read.table("wheat_depth.txt",header = F,sep = "\t")
head(data)
names(data) <- c('chr','start','stop','counts','atcg','length','depth')
head(data)
data$chr <- factor(data$chr,levels = c("chr1A","chr1B","chr1D","chr2A","chr2B","chr2D","chr3A","chr3B","chr3D","chr4A","chr4B","chr4D","chr5A","chr5B","chr5D","chr6A","chr6B","chr6D","chr7A","chr7B","chr7D","chrUn"))
#data$density <- c(log10(data$counts+1))
input <- read.table(" wheat_input_depth.txt",header = F,sep = "\t")
names(input) <- c('chr','start','stop','counts','atcg','length','depth')
data$nomalized <- data$counts/input$counts
head(data)
max(data$nomalized)
p1 <- ggplot(data,aes(start/1000000,nomalized,fill=chr),geom_area()+
    scale_y_continuous(breaks=c(0,5))+ylim(0,5)+
    alpha(0.5),legend(alpha(0.5),ncol = 1,cex = 0.6)+
    legend(x, y = NULL, fill = NULL,ncol = 1)+
    theme(axis.text.y.right = c("1A","1B","1D","2A","2B","2D","3A","3B","3D","4A","4B","4D","5A","5B","5D","6A","6B","6D","7A","7B","7D","Un"))+
    scale_y_discrete(labels(c("1A","1B","1D","2A","2B","2D","3A","3B","3D","4A","4B","4D","5A","5B","5D","6A","6B","6D","7A","7B","7D","Un"))))+
    facet_grid(chr~.,as.table = TRUE)+
    geom_bar(stat = "identity",width = .5)+
    scale_fill_manual(values =rainbow(22))+
    ylab("reads density")+xlab("chromosome position (Mb)")+labs(title = "Nomalized Reads Density in Chromosomes of wheat")
p1
ggsave(filename = "p1.pdf",height=10,width=6.18,plot=p1,dpi=600)

p2 
ggsave(filename = "p2.pdf",height=10,width=6.18,plot=p2,dpi=600)
               

效果图

未调整 ，大概就是这样