有关步骤可以参考之前的
2020-04-23 绘制Circos图之二 测序深度https://www.jianshu.com/p/ce76f043b45f
首先需要生成一个bed文件的滑窗
RNA转录组测序完成后，会进行序列比对，得到的结果是BAM格式的文件。
首先将生成的bam文件进行排序，
samtools sort -@ 10 -o ./sort.bam/19-1.sorted.bam ./bam/19-1.bam ###将bam文件排序
查看bam文件,格式
samtools view -H /home/qinghai/my_first_circos/sort.bam/19-1.sorted.bam
hg19.genome数据的分布必须与下图一致
可以网上ucsc下载hg19.genome
sudo apt install mysql-client-core-8.0#
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from hg19.chromInfo" > hg19.genome
结果如下，和bam文件序列不一致，别的物种也可以使用
http://genome.ucsc.edu/cgi-bin/hgGateway hgsid=1053922007_v51fhzoNz3YBAcyy1CIAVQtaeqM9
图片2.png
人用hg19，鼠mm10，斑马鱼danRer11
如斑马鱼：mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from danRer11
.chromInfo" > danRer11.genome
按照bam文件格式序列进行修改，如果哪个染色体不需要，直接改为0，在这里我们只需要线粒体，其他都改为0
bedtools makewindows 生成滑窗区域文本，可以理解-s -w 的区别，不清楚的可以直接输入
bedtools makewindows获取帮助信息。。。最好每次新生成一个划窗。
sudo apt install mysql-client-core-8.0#安装
bedtools makewindows -g hg19.genome -w 10000000 #间隔1000000
chr1 0 10000000
chr1 10000000 20000000
chr1 20000000 30000000
chr1 30000000 40000000
$ bedtools makewindows -g hg19.genome -w 1000000 -s 500000
chr1 0 1000000
chr1 500000 1500000
chr1 1000000 2000000
chr1 1500000 2500000
chr1 2000000 3000000
#bedtools makewindows -g hg19.genome -w 10000000 -s 1000000 > windows.bed
#or bedtools makewindows -g hg19.genome -w 100 -s 50 > windows.bed ##线粒体
bedtools makewindows -g hg19.genome -w 100 > windows.bed
> ### 注意需要删掉genome文件第一行, 要不然会一直报错
chrom size
chr4 78093715
chr7 74282399
chr5 72500376
chr3 62628489
chr6 60270059
chr2 59640629
...
结果如下，三列
bedtools makewindows用来自动生成划窗区间。-g genome.txt是要划分的基因组，格式为两列：染色体、染色体长度；-w 10000000为窗口大小为10M；-s 1000000为步长为1M，即窗口在染色体上每次向右平移1M的距离；windows.bed为输出文件，格式为三列：染色体、区间开始位点、区间结束位点。
自己编写hg19.txt或genome容易出错，需要通过mysql下载hg19.genome，然后根据需要进行修改才行。
报错。Error: chr4 has length equal to 0. Each chromosome must have non-zero length. Exiting.图片1.png
只需要将其他的染色体 改为 1就行，如图图片9.png
线粒体测序深度结果运行, 每隔多少bp进行深度测序，得到的数据可直接用于下游R分析绘图。
bedtools coverage -mean -sorted -g hg19.genome -a windows.bed -b /home/qinghai/my_first_circos/sort.bam/19-1.sorted.bam> 19-1.depth.txt
生成结果如下：图片12.png
这个结果还无法导入到mitochondrion.conf文件中，需要修改为下面的图。可通过excel编辑
karyotype.tair10.txt 中的信息 为，需要与之对应，因此第一列改为NC_027757.2图片12.png
chr - NC_027757.2 human_mito 0 16596 chr1该部分是用于circos 绘图
如果需要区分正负链
bedtools makewindows -g danio.genome -w 10 -i winnum > windows_MTO1.bed
需要加上-i winnum
-i winnum" - use the window number as the ID (e.g. 1,2,3,4...).
图片12.png
再通过excel编辑，在后面加上两列，特别是第六列改为+
再通过excel编辑，在后面加上两列，特别是第六列改为+
图片13.png
bedtools coverage -mean -sorted -g danio.genome -a windows_MTO1.bed -b ./sort_bam/Mmto1c1.sorted.bam -s> Mmto1c1.plus_depth.txt
###正链深度表
bedtools coverage -mean -sorted -g danio.genome -a windows_MTO1.bed -b ./sort_bam/Mmto1c1.sorted.bam -S> Mmto1c1.minus_depth.txt
###负链深度表
EXCEL中删除4，5，6列，并加上标题
用samtools对这个区间内的每个点做深度计算，代码如下：
samtools depth 19-1.sorted.bam -b windows.bed >test.stat
结果如下：结果如下所示（其中第一列为染色体，第二列为position，第三列即为这个position的深度）：
chrM 1 5708
chrM 2 5892
chrM 3 6009
chrM 4 6376
chrM 5 6516
chrM 6 6914
chrM 7 6977
。。。。。。。。。
chrM 16562 1873
chrM 16563 1851
chrM 16564 1803
chrM 16565 1744
chrM 16566 1659
chrM 16567 1620
chrM 16568 1568
chrM 16569 1473
》# 线粒体基因16569碱基，每个位点都进了深度测序统计。
这个是每个位点进行深度测序分析，目前不太会对间隔bp进行分析。
用shell对区间内所有的测序深度的depth做统计，代码如下：
$ cut -f3 test.stat |sort|uniq -c >test.stat.count
得到结果如下，其中第二列为depth，第一列为这个depth的频数，如第2行，即表示有1001个点的深度为2.
1 100
2 1001
1 1003
1 1004
1 1005
1 1007
接下来用R做这些深度分布图，代码如下：
》 首先用excel打开test.stat. 去除不需要的列数
1 5708
2 5892
3 6009
4 6376
5 6516
6 6914
7 6977
8 7135
9 7246
.....
也可以选择某一段区域进行作图，比如5062-8000区域
5062 3628
5063 3631
5064 3708
5065 3707
5066 3728
5067 3698
5068 3739
5069 3733
5070 3722
......... ##另存为test.stat2.csv
代码如下：
count<- read.csv('test.stat2.csv')
head(count)
names(count)<-c('X5062','X3628')###列的名字
ggplot(count,aes(x=X5062,y=X3628))+geom_point(col='red')
线图
ggplot(count,aes(x=v3,y=v4))+geom_line()
点线图
ggplot(count,aes(x=v3,y=v4))+geom_line()+geom_point() #默认颜色为黑色
基本格式：ggplot(data,aes(x,y))+geom_xx()+annotate()+labs()
data读取的excel数据，aes是美化的意思数据的列名，geom几何，annotate是注释，lab是标签
双折线图绘制
> count<- read.csv('MCK.depth.csv')
> head(count)
V1 V2
V1 V2 V3
1 4020 223.1 194.9
2 4030 224.6 195.8
3 4040 226.3 200.1
4 4050 218.1 195.2
5 4060 214.3 197.2
6 4070 211.3 196.2
需要对这一数据进行融合操作, 按照位置对数据进行融合之后，就得到了我们绘图所要的数据：
> library(reshape2)
> mydata<-melt(count,id="V1")
> head(mydata)
V1 variable value
1 4020 V2 223.1
2 4030 V2 224.6
3 4040 V2 226.3
4 4050 V2 218.1
5 4060 V2 214.3
6 4070 V2 211.3
融合之后的数据为3列，这里的三个列名是自动生成的，postion将默认为横轴title；variable将在绘图过程中默认为图例title，value将默认为纵轴title;如果需要对其进行改动，可在此处更改列名：
> colnames(mydata) <- c("position","sample","value")
> head(mydata)
position sample value
1 4020 V2 223.1
2 4030 V2 224.6
3 4040 V2 226.3
4 4050 V2 218.1
5 4060 V2 214.3
6 4070 V2 211.3
之后开始作图：
> ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+ geom_point()+ geom_line()
得到结果
进一步显示想要部分的区间的图。
ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+ geom_point()+ geom_line() + scale_x_continuous(limits = c(0000,10000),breaks = seq(5000,10000,500))
但是呢，这里存在一个问题，就是12S和16S rRNA量太多，导致其他部分量显示太少，就算只展示一部分区间的结果也不会更突出。因此在前期excel处理的时候可以把12S和16S数据先删除，在做5000-10000区间的数据图。
ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))
+ geom_point()+ geom_line()+labs(x="营业厅名称",y="评分",title="员工形象管理评分情况",fill="")
+theme(plot.title = element_text(hjust = 0.5))
ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))
+ geom_point()+ geom_line()+labs(x="营业厅名称",y="评分",title="员工形象管理评分情况",fill="")
+ theme(panel.grid.major=element_line(colour=NA), panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA), panel.grid.minor = element_blank(), ##以上theme中代码用于去除网格线且保留坐标轴边框
text = element_text(family = "STXihei"), #设置中文字体的显示
legend.position = c(.075,.915), #更改图例的位置，放至图内部的左上角
legend.box.background = element_rect(color="black")) #为图例田间边框线
还可以更改横坐标值
ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+ geom_point()+ geom_line() +labs(x="营业厅名称",y="评分",title="员工形象管理评分情况",fill="")+ theme(panel.grid.major=element_line(colour=NA),panel.background = element_rect(fill = "transparent",colour = NA),plot.background = element_rect(fill = "transparent",colour = NA), panel.grid.minor = element_blank(),text = element_text(family = "STXihei"),legend.position = c(.075,.915),legend.box.background = element_rect(color="black")) + scale_x_continuous(limits = c(0000,10000),breaks = seq(5000,10000,500))
其中V3为ck，V2为hom，前面可以改的，为了方便操作，就不改了
其实，用下面这一个作图的指令就够了
ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+ geom_line()+ scale_x_continuous(limits = c(4000,12000),breaks = seq(4000,12000,1000))
与线粒体序列对应着看，我们可以知道cox1 cox3 测序深度变低了
